diff --git a/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp b/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp index c18b720f28..bc66c7b026 100644 --- a/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp +++ b/Modules/DiffusionImaging/DicomImport/mitkDicomDiffusionImageHeaderReader.cpp @@ -1,345 +1,370 @@ /*========================================================================= 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::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 = 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::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 = m_VolumeReader->GetMetaDataDictionaryArray(); // load in all public tags 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 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], // 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 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 d925d42639..77adc61bfd 100644 --- a/Modules/DiffusionImaging/DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp +++ b/Modules/DiffusionImaging/DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp @@ -1,167 +1,192 @@ /*========================================================================= 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" #if GDCM_MAJOR_VERSION >= 2 #define DGDCM2 #endif #ifndef DGDCM2 #include "gdcm.h" // relevant Siemens private tags // relevant GE private tags const gdcm::DictEntry GEDictBValue( 0x0043, 0x1039, "IS", "1", "B Value of diffusion weighting" ); const gdcm::DictEntry GEDictXGradient( 0x0019, 0x10bb, "DS", "1", "X component of gradient direction" ); const gdcm::DictEntry GEDictYGradient( 0x0019, 0x10bc, "DS", "1", "Y component of gradient direction" ); const gdcm::DictEntry GEDictZGradient( 0x0019, 0x10bd, "DS", "1", "Z component of gradient direction" ); #else #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" ); #endif 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 = m_VolumeReader->GetMetaDataDictionaryArray(); #ifndef DGDCM2 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); #endif 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; 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; } 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], // vendor,SliceMosaic); diff --git a/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp b/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp index 44085c241a..445f1de484 100644 --- a/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp +++ b/Modules/DiffusionImaging/DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp @@ -1,667 +1,692 @@ /*========================================================================= 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 #if GDCM_MAJOR_VERSION >= 2 #define DGDCM2 #endif #ifndef DGDCM2 #include "gdcm.h" #else #include "gdcmFile.h" #include "gdcmImageReader.h" #include "gdcmDictEntry.h" #include "gdcmDicts.h" #include "gdcmTag.h" #endif 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; std::string::size_type pos = -1; while(nItems != 3) { 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 inputDict = m_VolumeReader->GetMetaDataDictionaryArray(); #ifndef DGDCM2 // relevant Siemens private tags gdcm::DictEntry SiemensDictBValue( 0x0019, 0x100c, "IS", "1", "B Value of diffusion weighting" ); gdcm::DictEntry SiemensDictDiffusionDirection( 0x0019, 0x100e, "FD", "3", "Diffusion Gradient Direction" ); gdcm::DictEntry SiemensDictDiffusionMatrix( 0x0019, 0x1027, "FD", "6", "Diffusion Matrix" ); gdcm::DictEntry SiemensDictShadowInfo( 0x0029, 0x1010, "OB", "1", "Siemens DWI Info" ); //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.GetName())==0) gdcm::Global::GetDicts()->GetDefaultPubDict()->AddEntry(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); #else // 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); #endif 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 ) { #ifndef DGDCM2 gdcm::File *header0 = new gdcm::File; gdcm::BinEntry* binEntry; header0->SetMaxSizeLoadEntry(65536); header0->SetFileName( m_DicomFilenames[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(); } #else 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)) { tag = std::string(ref.GetByteValue()->GetPointer(),ref.GetByteValue()->GetLength()); } } #endif // 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; #ifndef DGDCM2 success = itk::ExposeMetaData ( *(*inputDict)[0], "0019|100c", tag ); this->m_Output->bValue = atof( tag.c_str() ); #else 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; } } #endif tag.clear(); if(success) { if(this->m_Output->bValue == 0) { MITK_INFO << "BV: 0 (Baseline image)"; continue; } #ifndef DGDCM2 success = itk::ExposeMetaData ( *(*inputDict)[k], "0019|100e", tag); #else 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; } } #endif 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 { 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], // 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 // = 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 // { // 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); // 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); // vect3d.normalize(); // 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); // continue; // } // 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); // } // else // { // vect3d[0] = valueArray[0]; // vect3d[1] = valueArray[1]; // vect3d[2] = valueArray[2]; // DiffusionVectorsOrig.push_back(vect3d); // vect3d.normalize(); // 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], // vendor,SliceMosaic); // // // set m_Output // // } // //} // diff --git a/Modules/DiffusionImaging/Reconstruction/QuadProg.cpp b/Modules/DiffusionImaging/Reconstruction/QuadProg.cpp index 5831a457d0..3c1f715a0e 100644 --- a/Modules/DiffusionImaging/Reconstruction/QuadProg.cpp +++ b/Modules/DiffusionImaging/Reconstruction/QuadProg.cpp @@ -1,824 +1,839 @@ /* Author: Luca Di Gaspero DIEGM - University of Udine, Italy l.digaspero@uniud.it http://www.diegm.uniud.it/digaspero/ LICENSE This file is part of QuadProg++: a C++ library implementing the algorithm of Goldfarb and Idnani for the solution of a (convex) Quadratic Programming problem by means of an active-set dual method. Copyright (C) 2007-2009 Luca Di Gaspero. Copyright (C) 2009 Eric Moyer. QuadProg++ is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. QuadProg++ is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with QuadProg++. If not, see . */ #include #include #include #include #include #include #include "QuadProg.h" //#define TRACE_SOLVER #ifdef TRACE_SOLVER #undef TRACE_SOLVER #endif namespace QuadProgPP{ // Utility functions for updating some data needed by the solution method void compute_d(Vector& d, const Matrix& J, const Vector& np); void update_z(Vector& z, const Matrix& J, const Vector& d, int iq); void update_r(const Matrix& R, Vector& r, const Vector& d, int iq); bool add_constraint(Matrix& R, Matrix& J, Vector& d, int& iq, double& rnorm); void delete_constraint(Matrix& R, Matrix& J, Vector& A, Vector& u, int n, int p, int& iq, int l); // Utility functions for computing the Cholesky decomposition and solving // linear systems void cholesky_decomposition(Matrix& A); void cholesky_solve(const Matrix& L, Vector& x, const Vector& b); void forward_elimination(const Matrix& L, Vector& y, const Vector& b); void backward_elimination(const Matrix& U, Vector& x, const Vector& y); // Utility functions for computing the scalar product and the euclidean // distance between two numbers double scalar_product(const Vector& x, const Vector& y); double distance(double a, double b); // Utility functions for printing vectors and matrices void print_matrix(const char* name, const Matrix& A, int n = -1, int m = -1); template void print_vector(const char* name, const Vector& v, int n = -1); // The Solving function, implementing the Goldfarb-Idnani method double QuadProg::solve_quadprog(Matrix& G, Vector& g0, const Matrix& CE, const Vector& ce0, const Matrix& CI, const Vector& ci0, Vector& x) { std::ostringstream msg; + + std::locale C("C"); + std::locale originalLocale = msg.getloc(); + msg.imbue(C); + { //Ensure that the dimensions of the matrices and vectors can be //safely converted from unsigned int into to int without overflow. unsigned mx = std::numeric_limits::max(); if(G.ncols() >= mx || G.nrows() >= mx || CE.nrows() >= mx || CE.ncols() >= mx || CI.nrows() >= mx || CI.ncols() >= mx || ci0.size() >= mx || ce0.size() >= mx || g0.size() >= mx){ msg << "The dimensions of one of the input matrices or vectors were " << "too large." << std::endl << "The maximum allowable size for inputs to solve_quadprog is:" << mx << std::endl; throw std::logic_error(msg.str()); } } int n = G.ncols(), p = CE.ncols(), m = CI.ncols(); if ((int)G.nrows() != n) { msg << "The matrix G is not a square matrix (" << G.nrows() << " x " << G.ncols() << ")"; throw std::logic_error(msg.str()); } if ((int)CE.nrows() != n) { msg << "The matrix CE is incompatible (incorrect number of rows " << CE.nrows() << " , expecting " << n << ")"; throw std::logic_error(msg.str()); } if ((int)ce0.size() != p) { msg << "The vector ce0 is incompatible (incorrect dimension " << ce0.size() << ", expecting " << p << ")"; throw std::logic_error(msg.str()); } if ((int)CI.nrows() != n) { msg << "The matrix CI is incompatible (incorrect number of rows " << CI.nrows() << " , expecting " << n << ")"; throw std::logic_error(msg.str()); } if ((int)ci0.size() != m) { msg << "The vector ci0 is incompatible (incorrect dimension " << ci0.size() << ", expecting " << m << ")"; throw std::logic_error(msg.str()); } x.resize(n); register int i, j, k, l; /* indices */ int ip; // this is the index of the constraint to be added to the active set Matrix R(n, n), J(n, n); Vector s(m + p), z(n), r(m + p), d(n), np(n), u(m + p), x_old(n), u_old(m + p); double f_value, psi, c1, c2, sum, ss, R_norm; double inf; if (std::numeric_limits::has_infinity) inf = std::numeric_limits::infinity(); else inf = 1.0E300; double t, t1, t2; /* t is the step lenght, which is the minimum of the partial step length t1 * and the full step length t2 */ Vector A(m + p), A_old(m + p), iai(m + p); int q, iq, iter = 0; Vector iaexcl(m + p); /* p is the number of equality constraints */ /* m is the number of inequality constraints */ q = 0; /* size of the active set A (containing the indices of the active constraints) */ #ifdef TRACE_SOLVER std::cout << std::endl << "Starting solve_quadprog" << std::endl; print_matrix("G", G); print_vector("g0", g0); print_matrix("CE", CE); print_vector("ce0", ce0); print_matrix("CI", CI); print_vector("ci0", ci0); #endif /* * Preprocessing phase */ /* compute the trace of the original matrix G */ c1 = 0.0; for (i = 0; i < n; i++) { c1 += G[i][i]; } /* decompose the matrix G in the form L^T L */ cholesky_decomposition(G); #ifdef TRACE_SOLVER print_matrix("G", G); #endif /* initialize the matrix R */ for (i = 0; i < n; i++) { d[i] = 0.0; for (j = 0; j < n; j++) R[i][j] = 0.0; } R_norm = 1.0; /* this variable will hold the norm of the matrix R */ /* compute the inverse of the factorized matrix G^-1, this is the initial value for H */ c2 = 0.0; for (i = 0; i < n; i++) { d[i] = 1.0; forward_elimination(G, z, d); for (j = 0; j < n; j++) J[i][j] = z[j]; c2 += z[i]; d[i] = 0.0; } #ifdef TRACE_SOLVER print_matrix("J", J); #endif /* c1 * c2 is an estimate for cond(G) */ /* * Find the unconstrained minimizer of the quadratic form 0.5 * x G x + g0 x * this is a feasible point in the dual space * x = G^-1 * g0 */ cholesky_solve(G, x, g0); for (i = 0; i < n; i++) x[i] = -x[i]; /* and compute the current solution value */ f_value = 0.5 * scalar_product(g0, x); #ifdef TRACE_SOLVER std::cout << "Unconstrained solution: " << f_value << std::endl; print_vector("x", x); #endif /* Add equality constraints to the working set A */ iq = 0; for (i = 0; i < p; i++) { for (j = 0; j < n; j++) np[j] = CE[j][i]; compute_d(d, J, np); update_z(z, J, d, iq); update_r(R, r, d, iq); #ifdef TRACE_SOLVER print_matrix("R", R, n, iq); print_vector("z", z); print_vector("r", r, iq); print_vector("d", d); #endif /* compute full step length t2: i.e., the minimum step in primal space s.t. the contraint becomes feasible */ t2 = 0.0; if (fabs(scalar_product(z, z)) > std::numeric_limits::epsilon()) // i.e. z != 0 t2 = (-scalar_product(np, x) - ce0[i]) / scalar_product(z, np); /* set x = x + t2 * z */ for (k = 0; k < n; k++) x[k] += t2 * z[k]; /* set u = u+ */ u[iq] = t2; for (k = 0; k < iq; k++) u[k] -= t2 * r[k]; /* compute the new solution value */ f_value += 0.5 * (t2 * t2) * scalar_product(z, np); A[i] = -i - 1; if (!add_constraint(R, J, d, iq, R_norm)) { // Equality constraints are linearly dependent throw std::runtime_error("Constraints are linearly dependent"); return f_value; } } /* set iai = K \ A */ for (i = 0; i < m; i++) iai[i] = i; l1: iter++; #ifdef TRACE_SOLVER print_vector("x", x); #endif /* step 1: choose a violated constraint */ for (i = p; i < iq; i++) { ip = A[i]; iai[ip] = -1; } /* compute s[x] = ci^T * x + ci0 for all elements of K \ A */ ss = 0.0; psi = 0.0; /* this value will contain the sum of all infeasibilities */ ip = 0; /* ip will be the index of the chosen violated constraint */ for (i = 0; i < m; i++) { iaexcl[i] = true; sum = 0.0; for (j = 0; j < n; j++) sum += CI[j][i] * x[j]; sum += ci0[i]; s[i] = sum; psi += std::min(0.0, sum); } #ifdef TRACE_SOLVER print_vector("s", s, m); #endif if (fabs(psi) <= m * std::numeric_limits::epsilon() * c1 * c2* 100.0) { /* numerically there are not infeasibilities anymore */ q = iq; return f_value; } /* save old values for u and A */ for (i = 0; i < iq; i++) { u_old[i] = u[i]; A_old[i] = A[i]; } /* and for x */ for (i = 0; i < n; i++) x_old[i] = x[i]; l2: /* Step 2: check for feasibility and determine a new S-pair */ for (i = 0; i < m; i++) { if (s[i] < ss && iai[i] != -1 && iaexcl[i]) { ss = s[i]; ip = i; } } if (ss >= 0.0) { q = iq; return f_value; } /* set np = n[ip] */ for (i = 0; i < n; i++) np[i] = CI[i][ip]; /* set u = [u 0]^T */ u[iq] = 0.0; /* add ip to the active set A */ A[iq] = ip; #ifdef TRACE_SOLVER std::cout << "Trying with constraint " << ip << std::endl; print_vector("np", np); #endif l2a:/* Step 2a: determine step direction */ /* compute z = H np: the step direction in the primal space (through J, see the paper) */ compute_d(d, J, np); update_z(z, J, d, iq); /* compute N* np (if q > 0): the negative of the step direction in the dual space */ update_r(R, r, d, iq); #ifdef TRACE_SOLVER std::cout << "Step direction z" << std::endl; print_vector("z", z); print_vector("r", r, iq + 1); print_vector("u", u, iq + 1); print_vector("d", d); print_vector("A", A, iq + 1); #endif /* Step 2b: compute step length */ l = 0; /* Compute t1: partial step length (maximum step in dual space without violating dual feasibility */ t1 = inf; /* +inf */ /* find the index l s.t. it reaches the minimum of u+[x] / r */ for (k = p; k < iq; k++) { if (r[k] > 0.0) { if (u[k] / r[k] < t1) { t1 = u[k] / r[k]; l = A[k]; } } } /* Compute t2: full step length (minimum step in primal space such that the constraint ip becomes feasible */ if (fabs(scalar_product(z, z)) > std::numeric_limits::epsilon()) // i.e. z != 0 t2 = -s[ip] / scalar_product(z, np); else t2 = inf; /* +inf */ /* the step is chosen as the minimum of t1 and t2 */ t = std::min(t1, t2); #ifdef TRACE_SOLVER std::cout << "Step sizes: " << t << " (t1 = " << t1 << ", t2 = " << t2 << ") "; #endif /* Step 2c: determine new S-pair and take step: */ /* case (i): no step in primal or dual space */ if (t >= inf) { /* QPP is infeasible */ // FIXME: unbounded to raise q = iq; return inf; } /* case (ii): step in dual space */ if (t2 >= inf) { /* set u = u + t * [-r 1] and drop constraint l from the active set A */ for (k = 0; k < iq; k++) u[k] -= t * r[k]; u[iq] += t; iai[l] = l; delete_constraint(R, J, A, u, n, p, iq, l); #ifdef TRACE_SOLVER std::cout << " in dual space: " << f_value << std::endl; print_vector("x", x); print_vector("z", z); print_vector("A", A, iq + 1); #endif goto l2a; } /* case (iii): step in primal and dual space */ /* set x = x + t * z */ for (k = 0; k < n; k++) x[k] += t * z[k]; /* update the solution value */ f_value += t * scalar_product(z, np) * (0.5 * t + u[iq]); /* u = u + t * [-r 1] */ for (k = 0; k < iq; k++) u[k] -= t * r[k]; u[iq] += t; #ifdef TRACE_SOLVER std::cout << " in both spaces: " << f_value << std::endl; print_vector("x", x); print_vector("u", u, iq + 1); print_vector("r", r, iq + 1); print_vector("A", A, iq + 1); #endif if (fabs(t - t2) < std::numeric_limits::epsilon()) { #ifdef TRACE_SOLVER std::cout << "Full step has taken " << t << std::endl; print_vector("x", x); #endif /* full step has taken */ /* add constraint ip to the active set*/ if (!add_constraint(R, J, d, iq, R_norm)) { iaexcl[ip] = false; delete_constraint(R, J, A, u, n, p, iq, ip); #ifdef TRACE_SOLVER print_matrix("R", R); print_vector("A", A, iq); print_vector("iai", iai); #endif for (i = 0; i < m; i++) iai[i] = i; for (i = p; i < iq; i++) { A[i] = A_old[i]; u[i] = u_old[i]; iai[A[i]] = -1; } for (i = 0; i < n; i++) x[i] = x_old[i]; goto l2; /* go to step 2 */ } else iai[ip] = -1; #ifdef TRACE_SOLVER print_matrix("R", R); print_vector("A", A, iq); print_vector("iai", iai); #endif goto l1; } /* a patial step has taken */ #ifdef TRACE_SOLVER std::cout << "Partial step has taken " << t << std::endl; print_vector("x", x); #endif /* drop constraint l */ iai[l] = l; delete_constraint(R, J, A, u, n, p, iq, l); #ifdef TRACE_SOLVER print_matrix("R", R); print_vector("A", A, iq); #endif /* update s[ip] = CI * x + ci0 */ sum = 0.0; for (k = 0; k < n; k++) sum += CI[k][ip] * x[k]; s[ip] = sum + ci0[ip]; #ifdef TRACE_SOLVER print_vector("s", s, m); #endif goto l2a; + + msg.imbue( originalLocale ); } inline void compute_d(Vector& d, const Matrix& J, const Vector& np) { register int i, j, n = d.size(); register double sum; /* compute d = H^T * np */ for (i = 0; i < n; i++) { sum = 0.0; for (j = 0; j < n; j++) sum += J[j][i] * np[j]; d[i] = sum; } } inline void update_z(Vector& z, const Matrix& J, const Vector& d, int iq) { register int i, j, n = z.size(); /* setting of z = H * d */ for (i = 0; i < n; i++) { z[i] = 0.0; for (j = iq; j < n; j++) z[i] += J[i][j] * d[j]; } } inline void update_r(const Matrix& R, Vector& r, const Vector& d, int iq) { register int i, j;/*, n = d.size();*/ register double sum; /* setting of r = R^-1 d */ for (i = iq - 1; i >= 0; i--) { sum = 0.0; for (j = i + 1; j < iq; j++) sum += R[i][j] * r[j]; r[i] = (d[i] - sum) / R[i][i]; } } bool add_constraint(Matrix& R, Matrix& J, Vector& d, int& iq, double& R_norm) { int n = d.size(); #ifdef TRACE_SOLVER std::cout << "Add constraint " << iq << '/'; #endif register int i, j, k; double cc, ss, h, t1, t2, xny; /* we have to find the Givens rotation which will reduce the element d[j] to zero. if it is already zero we don't have to do anything, except of decreasing j */ for (j = n - 1; j >= iq + 1; j--) { /* The Givens rotation is done with the matrix (cc cs, cs -cc). If cc is one, then element (j) of d is zero compared with element (j - 1). Hence we don't have to do anything. If cc is zero, then we just have to switch column (j) and column (j - 1) of J. Since we only switch columns in J, we have to be careful how we update d depending on the sign of gs. Otherwise we have to apply the Givens rotation to these columns. The i - 1 element of d has to be updated to h. */ cc = d[j - 1]; ss = d[j]; h = distance(cc, ss); if (fabs(h) < std::numeric_limits::epsilon()) // h == 0 continue; d[j] = 0.0; ss = ss / h; cc = cc / h; if (cc < 0.0) { cc = -cc; ss = -ss; d[j - 1] = -h; } else d[j - 1] = h; xny = ss / (1.0 + cc); for (k = 0; k < n; k++) { t1 = J[k][j - 1]; t2 = J[k][j]; J[k][j - 1] = t1 * cc + t2 * ss; J[k][j] = xny * (t1 + J[k][j - 1]) - t2; } } /* update the number of constraints added*/ iq++; /* To update R we have to put the iq components of the d vector into column iq - 1 of R */ for (i = 0; i < iq; i++) R[i][iq - 1] = d[i]; #ifdef TRACE_SOLVER std::cout << iq << std::endl; print_matrix("R", R, iq, iq); print_matrix("J", J); print_vector("d", d, iq); #endif if (fabs(d[iq - 1]) <= std::numeric_limits::epsilon() * R_norm) { // problem degenerate return false; } R_norm = std::max(R_norm, fabs(d[iq - 1])); return true; } void delete_constraint(Matrix& R, Matrix& J, Vector& A, Vector& u, int n, int p, int& iq, int l) { #ifdef TRACE_SOLVER std::cout << "Delete constraint " << l << ' ' << iq; #endif register int i, j, k, qq = -1; // just to prevent warnings from smart compilers double cc, ss, h, xny, t1, t2; /* Find the index qq for active constraint l to be removed */ for (i = p; i < iq; i++) if (A[i] == l) { qq = i; break; } /* remove the constraint from the active set and the duals */ for (i = qq; i < iq - 1; i++) { A[i] = A[i + 1]; u[i] = u[i + 1]; for (j = 0; j < n; j++) R[j][i] = R[j][i + 1]; } A[iq - 1] = A[iq]; u[iq - 1] = u[iq]; A[iq] = 0; u[iq] = 0.0; for (j = 0; j < iq; j++) R[j][iq - 1] = 0.0; /* constraint has been fully removed */ iq--; #ifdef TRACE_SOLVER std::cout << '/' << iq << std::endl; #endif if (iq == 0) return; for (j = qq; j < iq; j++) { cc = R[j][j]; ss = R[j + 1][j]; h = distance(cc, ss); if (fabs(h) < std::numeric_limits::epsilon()) // h == 0 continue; cc = cc / h; ss = ss / h; R[j + 1][j] = 0.0; if (cc < 0.0) { R[j][j] = -h; cc = -cc; ss = -ss; } else R[j][j] = h; xny = ss / (1.0 + cc); for (k = j + 1; k < iq; k++) { t1 = R[j][k]; t2 = R[j + 1][k]; R[j][k] = t1 * cc + t2 * ss; R[j + 1][k] = xny * (t1 + R[j][k]) - t2; } for (k = 0; k < n; k++) { t1 = J[k][j]; t2 = J[k][j + 1]; J[k][j] = t1 * cc + t2 * ss; J[k][j + 1] = xny * (J[k][j] + t1) - t2; } } } inline double distance(double a, double b) { register double a1, b1, t; a1 = fabs(a); b1 = fabs(b); if (a1 > b1) { t = (b1 / a1); return a1 * ::std::sqrt(1.0 + t * t); } else if (b1 > a1) { t = (a1 / b1); return b1 * ::std::sqrt(1.0 + t * t); } return a1 * ::std::sqrt(2.0); } inline double scalar_product(const Vector& x, const Vector& y) { register int i, n = x.size(); register double sum; sum = 0.0; for (i = 0; i < n; i++) sum += x[i] * y[i]; return sum; } void cholesky_decomposition(Matrix& A) { register int i, j, k, n = A.nrows(); register double sum; + std::locale C("C"); for (i = 0; i < n; i++) { for (j = i; j < n; j++) { sum = A[i][j]; for (k = i - 1; k >= 0; k--) sum -= A[i][k]*A[j][k]; if (i == j) { if (sum <= 0.0) { std::ostringstream os; + os.imbue(C); // raise error print_matrix("A", A); os << "Error in cholesky decomposition, sum: " << sum; throw std::logic_error(os.str()); exit(-1); } A[i][i] = ::std::sqrt(sum); } else A[j][i] = sum / A[i][i]; } for (k = i + 1; k < n; k++) A[i][k] = A[k][i]; } } void cholesky_solve(const Matrix& L, Vector& x, const Vector& b) { int n = L.nrows(); Vector y(n); /* Solve L * y = b */ forward_elimination(L, y, b); /* Solve L^T * x = y */ backward_elimination(L, x, y); } inline void forward_elimination(const Matrix& L, Vector& y, const Vector& b) { register int i, j, n = L.nrows(); y[0] = b[0] / L[0][0]; for (i = 1; i < n; i++) { y[i] = b[i]; for (j = 0; j < i; j++) y[i] -= L[i][j] * y[j]; y[i] = y[i] / L[i][i]; } } inline void backward_elimination(const Matrix& U, Vector& x, const Vector& y) { register int i, j, n = U.nrows(); x[n - 1] = y[n - 1] / U[n - 1][n - 1]; for (i = n - 2; i >= 0; i--) { x[i] = y[i]; for (j = i + 1; j < n; j++) x[i] -= U[i][j] * x[j]; x[i] = x[i] / U[i][i]; } } void print_matrix(const char* name, const Matrix& A, int n, int m) { std::ostringstream s; + std::locale C("C"); + s.imbue(C); + std::string t; if (n == -1) n = A.nrows(); if (m == -1) m = A.ncols(); s << name << ": " << std::endl; for (int i = 0; i < n; i++) { s << " "; for (int j = 0; j < m; j++) s << A[i][j] << ", "; s << std::endl; } t = s.str(); t = t.substr(0, t.size() - 3); // To remove the trailing space, comma and newline std::cout << t << std::endl; } template void print_vector(const char* name, const Vector& v, int n) { std::ostringstream s; + std::locale C("C"); + s.imbue(C); + std::string t; if (n == -1) n = v.size(); s << name << ": " << std::endl << " "; for (int i = 0; i < n; i++) { s << v[i] << ", "; } t = s.str(); t = t.substr(0, t.size() - 2); // To remove the trailing space and comma std::cout << t << std::endl; } } diff --git a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp index 511cd4a790..a91cdfb06b 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp @@ -1,884 +1,897 @@ /*========================================================================= Program: Insight Segmentation & Registration Toolkit Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.txx,v $ Language: C++ Date: $Date: 2006-07-19 15:11:41 $ Version: $Revision: 1.11 $ Copyright (c) Insight Software Consortium. All rights reserved. See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp #define __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp -#include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" -#include "itkImageRegionConstIterator.h" -#include "itkImageRegionConstIteratorWithIndex.h" -#include "itkImageRegionIterator.h" -#include "itkArray.h" -#include "vnl/vnl_vector.h" +#include +#include +#include +#include +#include +#include #include +#include +#include #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 #include #endif #endif #include "itkPointShell.h" namespace itk { #define QBALL_ANAL_RECON_PI 3.14159265358979323846 template< class T, class TG, class TO, int L, int NODF> AnalyticalDiffusionQballReconstructionImageFilter ::AnalyticalDiffusionQballReconstructionImageFilter() : m_GradientDirectionContainer(NULL), m_NumberOfGradientDirections(0), m_NumberOfBaselineImages(1), m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), m_BValue(1.0), m_Lambda(0.0), m_DirectionsDuplicated(false) { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> typename itk::AnalyticalDiffusionQballReconstructionImageFilter< TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType, NOrderL,NrOdfDirections>::OdfPixelType itk::AnalyticalDiffusionQballReconstructionImageFilter ::Normalize( OdfPixelType odf, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { case QBAR_STANDARD: { TOdfPixelType sum = 0; for(int i=0; i0) odf /= sum; return odf; break; } case QBAR_B_ZERO_B_VALUE: { for(int i=0; i0) odf /= sum; break; } case QBAR_NONNEG_SOLID_ANGLE: { break; } } return odf; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> vnl_vector itk::AnalyticalDiffusionQballReconstructionImageFilter ::PreNormalize( vnl_vector vec, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { case QBAR_STANDARD: { return vec; break; } case QBAR_B_ZERO_B_VALUE: { int n = vec.size(); for(int i=0; i= b0f) meas = b0f - 0.001; vec[i] = log(-log(meas/b0f)); } return vec; break; } } return vec; } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::BeforeThreadedGenerateData() { // If we have more than 2 inputs, then each input, except the first is a // gradient image. The number of gradient images must match the number of // gradient directions. //const unsigned int numberOfInputs = this->GetNumberOfInputs(); // There need to be at least 6 gradient directions to be able to compute the // tensor basis if( m_NumberOfGradientDirections < 6 ) { itkExceptionMacro( << "At least 6 gradient directions are required" ); } // Input must be an itk::VectorImage. std::string gradientImageClassName( this->ProcessObject::GetInput(0)->GetNameOfClass()); if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) { itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. " << "But its of type: " << gradientImageClassName ); } this->ComputeReconstructionMatrix(); m_BZeroImage = BZeroImageType::New(); typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_BZeroImage->Allocate(); m_ODFSumImage = BZeroImageType::New(); m_ODFSumImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_ODFSumImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_ODFSumImage->SetDirection( img->GetDirection() ); // Set the image direction m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); m_ODFSumImage->Allocate(); if(m_NormalizationMethod == QBAR_SOLID_ANGLE || m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) { m_Lambda = 0.0; } } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int ) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); oit2.GoToBegin(); ImageRegionIterator< BlaImage > oit3(m_ODFSumImage, outputRegionForThread); oit3.GoToBegin(); typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; typedef typename GradientImagesType::PixelType GradientVectorType; typename GradientImagesType::Pointer gradientImagePointer = NULL; // Would have liked a dynamic_cast here, but seems SGI doesn't like it // The enum will ensure that an inappropriate cast is not done gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); GradientIteratorType git(gradientImagePointer, outputRegionForThread ); git.GoToBegin(); // Compute the indicies of the baseline images and gradient images std::vector baselineind; // contains the indicies of // the baseline images std::vector gradientind; // contains the indicies of // the gradient images for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() <= 0.0) { baselineind.push_back(gdcit.Index()); } else { gradientind.push_back(gdcit.Index()); } } if( m_DirectionsDuplicated ) { int gradIndSize = gradientind.size(); for(int i=0; i::AccumulateType b0 = NumericTraits::Zero; // Average the baseline image pixels for(unsigned int i = 0; i < baselineind.size(); ++i) { b0 += b[baselineind[i]]; } b0 /= this->m_NumberOfBaselineImages; OdfPixelType odf(0.0); vnl_vector B(m_NumberOfGradientDirections); if( (b0 != 0) && (b0 >= m_Threshold) ) { for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { B[i] = static_cast(b[gradientind[i]]); } B = PreNormalize(B, b0); if(m_NormalizationMethod == QBAR_SOLID_ANGLE) { vnl_vector coeffs(m_NumberCoefficients); coeffs = ( (*m_CoeffReconstructionMatrix) * B ); coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block(); } else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE) { /** this would be the place to implement a non-negative * solver for quadratic programming problem: * min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones * (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints") * .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's */ QuadProgPP::Matrix& G(m_G); int lenb = B.size(); vnl_vector* s = new vnl_vector(lenb); for (int ii=0; ii g0_ = -1.0 * (*m_B_t) * (*s); QuadProgPP::Vector g0(g0_.data_block(),m_NumberCoefficients); try { QuadProgPP::QuadProg::solve_quadprog(G,g0,m_CE,m_ce0,m_CI,m_ci0,m_x); } catch(...) { m_x = 0; } vnl_vector coeffs(m_NumberCoefficients); for (int ii=0; ii void AnalyticalDiffusionQballReconstructionImageFilter ::tofile2(vnl_matrix *pA, std::string fname) { vnl_matrix A = (*pA); ofstream myfile; + std::locale C("C"); + std::locale originalLocale = myfile.getloc(); + myfile.imbue(C); + myfile.open (fname.c_str()); myfile << "A1=["; for(int i=0; i double AnalyticalDiffusionQballReconstructionImageFilter ::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::Cart2Sph(double x, double y, double z, double *cart) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); th = atan2(z,sqrt(x*x+y*y)); phi = atan2(y,x); th = -th+QBALL_ANAL_RECON_PI/2; phi = -phi+QBALL_ANAL_RECON_PI; cart[0] = phi; cart[1] = th; cart[2] = rad; } template< class T, class TG, class TO, int L, int NODF> double AnalyticalDiffusionQballReconstructionImageFilter ::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i double AnalyticalDiffusionQballReconstructionImageFilter ::spherical_harmonic(int m,int l,double theta,double phi, bool complexPart) { if( theta < 0 || theta > QBALL_ANAL_RECON_PI ) { std::cout << "bad range" << std::endl; return 0; } if( phi < 0 || phi > 2*QBALL_ANAL_RECON_PI ) { std::cout << "bad range" << std::endl; return 0; } double pml = 0; double fac1 = factorial(l+m); double fac2 = factorial(l-m); if( m<0 ) { #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 pml = ::boost::math::legendre_p(l, -m, cos(theta)); #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; #endif #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; #endif double mypow = pow(-1.0,-m); double myfac = (fac1/fac2); pml *= mypow*myfac; } else { #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 pml = ::boost::math::legendre_p(l, m, cos(theta)); #endif #endif } //std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl; double retval = sqrt(((2.0*(double)l+1.0)/(4.0*QBALL_ANAL_RECON_PI))*(fac2/fac1)) * pml; if( !complexPart ) { retval *= cos(m*phi); } else { retval *= sin(m*phi); } //std::cout << retval << std::endl; return retval; } template< class T, class TG, class TO, int L, int NODF> double AnalyticalDiffusionQballReconstructionImageFilter ::Yj(int m, int k, double theta, double phi) { if( -k <= m && m < 0 ) { return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,false); } if( m == 0 ) return spherical_harmonic(0,k,theta,phi,false); if( 0 < m && m <= k ) { return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true); } return 0; } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::ComputeReconstructionMatrix() { //for(int i=-6;i<7;i++) // std::cout << boost::math::legendre_p(6, i, 0.65657) << std::endl; if( m_NumberOfGradientDirections < 6 ) { itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" ); } { // check for duplicate diffusion gradients bool warning = false; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) { if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index()) { itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); warning = true; break; } } if (warning) break; } // handle acquisition schemes where only half of the spherical // shell is sampled by the gradient directions. In this case, // each gradient direction is duplicated in negative direction. vnl_vector centerMass(3); centerMass.fill(0.0); int count = 0; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { if(gdcit1.Value().one_norm() > 0.0) { centerMass += gdcit1.Value(); count ++; } } centerMass /= count; if(centerMass.two_norm() > 0.1) { m_DirectionsDuplicated = true; m_NumberOfGradientDirections *= 2; } } vnl_matrix *Q = new vnl_matrix(3, m_NumberOfGradientDirections); { int i = 0; for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { double x = gdcit.Value().get(0); double y = gdcit.Value().get(1); double z = gdcit.Value().get(2); double cart[3]; Cart2Sph(x,y,z,cart); (*Q)(0,i) = cart[0]; (*Q)(1,i) = cart[1]; (*Q)(2,i++) = cart[2]; } } if(m_DirectionsDuplicated) { for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { double x = gdcit.Value().get(0); double y = gdcit.Value().get(1); double z = gdcit.Value().get(2); double cart[3]; Cart2Sph(x,y,z,cart); (*Q)(0,i) = cart[0]; (*Q)(1,i) = cart[1]; (*Q)(2,i++) = cart[2]; } } } } int l = L; m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l; vnl_matrix* B = new vnl_matrix(m_NumberOfGradientDirections,m_NumberCoefficients); vnl_matrix* _L = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); _L->fill(0.0); vnl_matrix* LL = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); LL->fill(0.0); vnl_matrix* P = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); P->fill(0.0); vnl_matrix* Inv = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients); P->fill(0.0); vnl_vector* lj = new vnl_vector(m_NumberCoefficients); m_LP = new vnl_vector(m_NumberCoefficients); for(unsigned int i=0; i(B->transpose()); //tofile2(&m_B_t,"m_B_t"); vnl_matrix B_t_B = (*m_B_t) * (*B); //tofile2(&B_t_B,"B_t_B"); vnl_matrix lambdaLL(m_NumberCoefficients,m_NumberCoefficients); lambdaLL.update((*LL)); lambdaLL *= m_Lambda; //tofile2(&lambdaLL,"lLL"); vnl_matrix tmp( B_t_B + lambdaLL); vnl_matrix_inverse *pseudoInverse = new vnl_matrix_inverse( tmp ); (*Inv) = pseudoInverse->pinverse(); //tofile2(Inv,"Inv"); vnl_matrix temp((*Inv) * (*m_B_t)); double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); switch(m_NormalizationMethod) { case QBAR_ADC_ONLY: case QBAR_RAW_SIGNAL: break; case QBAR_STANDARD: case QBAR_B_ZERO_B_VALUE: case QBAR_B_ZERO: case QBAR_NONE: temp = (*P) * temp; break; case QBAR_SOLID_ANGLE: temp = fac1 * (*P) * (*_L) * temp; break; case QBAR_NONNEG_SOLID_ANGLE: m_G = QuadProgPP::Matrix(B_t_B.data_block(), B_t_B.rows(), B_t_B.cols()); m_CE = QuadProgPP::Matrix((double)0,m_NumberCoefficients,0); m_ce0 = QuadProgPP::Vector((double)0,0); m_ci0 = QuadProgPP::Vector(4*QBALL_ANAL_RECON_PI, NODF); m_x = QuadProgPP::Vector(m_NumberCoefficients); (*m_LP) *= fac1; break; } //tofile2(&temp,"A"); m_CoeffReconstructionMatrix = new vnl_matrix(m_NumberCoefficients,m_NumberOfGradientDirections); for(int i=0; iodfs later int NOdfDirections = NODF; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); m_SphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); vnl_matrix* sphericalHarmonicBasisMatrix2 = new vnl_matrix(NOdfDirections,m_NumberCoefficients); for(int i=0; i CI_t = (*sphericalHarmonicBasisMatrix2) * (*P) * (*_L); m_CI = QuadProgPP::Matrix(CI_t.transpose().data_block(), m_NumberCoefficients, NOdfDirections); } m_ReconstructionMatrix = new vnl_matrix(NOdfDirections,m_NumberOfGradientDirections); *m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix); } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { this->m_GradientDirectionContainer = gradientDirection; unsigned int numImages = gradientDirection->Size(); this->m_NumberOfBaselineImages = 0; for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin(); it != this->m_GradientDirectionContainer->End(); it++) { if(it.Value().one_norm() <= 0.0) { this->m_NumberOfBaselineImages++; } else // Normalize non-zero gradient directions { it.Value() = it.Value() / it.Value().two_norm(); } } this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; // ensure that the gradient image we received has as many components as // the number of gradient directions if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) { itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages << " directions specified but image has " << gradientImage->GetVectorLength() << " components."); } this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); } template< class T, class TG, class TO, int L, int NODF> void AnalyticalDiffusionQballReconstructionImageFilter ::PrintSelf(std::ostream& os, Indent indent) const { + std::locale C("C"); + std::locale originalLocale = os.getloc(); + os.imbue(C); + Superclass::PrintSelf(os,indent); os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl; if ( m_GradientDirectionContainer ) { os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl; } else { os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl; } os << indent << "NumberOfGradientDirections: " << m_NumberOfGradientDirections << std::endl; os << indent << "NumberOfBaselineImages: " << m_NumberOfBaselineImages << std::endl; os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl; os << indent << "BValue: " << m_BValue << std::endl; + + os.imbue( originalLocale ); } } #endif // __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx b/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx index 112091befb..5f1fe59525 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionQballReconstructionImageFilter.txx @@ -1,857 +1,863 @@ /*========================================================================= Program: Insight Segmentation & Registration Toolkit Language: C++ Date: $Date: 2006-07-19 15:11:41 $ Version: $Revision: 1.11 $ Copyright (c) Insight Software Consortium. All rights reserved. See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 __itkDiffusionQballReconstructionImageFilter_txx #define __itkDiffusionQballReconstructionImageFilter_txx #include "itkDiffusionQballReconstructionImageFilter.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionIterator.h" #include "itkArray.h" #include "vnl/vnl_vector.h" #include "itkPointShell.h" namespace itk { #define QBALL_RECON_PI 3.14159265358979323846 template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::DiffusionQballReconstructionImageFilter() : m_GradientDirectionContainer(NULL), m_NumberOfGradientDirections(0), m_NumberOfEquatorSamplingPoints(0), m_NumberOfBaselineImages(1), m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), m_BValue(1.0), m_GradientImageTypeEnumeration(Else), m_DirectionsDuplicated(false) { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::BeforeThreadedGenerateData() { // If we have more than 2 inputs, then each input, except the first is a // gradient image. The number of gradient images must match the number of // gradient directions. const unsigned int numberOfInputs = this->GetNumberOfInputs(); // There need to be at least 6 gradient directions to be able to compute the // tensor basis if( m_NumberOfGradientDirections < 6 ) { itkExceptionMacro( << "At least 6 gradient directions are required" ); } // If there is only 1 gradient image, it must be an itk::VectorImage. Otherwise // we must have a container of (numberOfInputs-1) itk::Image. Check to make sure if ( numberOfInputs == 1 && m_GradientImageTypeEnumeration != GradientIsInASingleImage ) { std::string gradientImageClassName( this->ProcessObject::GetInput(0)->GetNameOfClass()); if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) { itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. " << "But its of type: " << gradientImageClassName ); } } // Compute reconstruction matrix that is multiplied to the data-vector // each voxel in order to reconstruct the ODFs this->ComputeReconstructionMatrix(); // Allocate the b-zero image m_BZeroImage = BZeroImageType::New(); if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) { typename ReferenceImageType::Pointer img = static_cast< ReferenceImageType * >(this->ProcessObject::GetInput(0)); m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); } // The gradients are specified in a single multi-component image else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage ) { typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); } m_BZeroImage->Allocate(); } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> typename itk::DiffusionQballReconstructionImageFilter::OdfPixelType itk::DiffusionQballReconstructionImageFilter::Normalize( OdfPixelType odf, typename NumericTraits::AccumulateType b0 ) { switch( m_NormalizationMethod ) { // divide by sum to retreive a PDF case QBR_STANDARD: { odf.Normalize(); return odf; break; } // ADC style case QBR_B_ZERO_B_VALUE: { for(int i=0; i vnl_vector itk::DiffusionQballReconstructionImageFilter::PreNormalize( vnl_vector vec ) { switch( m_NormalizationMethod ) { // standard: no normalization before reconstruction case QBR_STANDARD: { return vec; break; } // log of signal case QBR_B_ZERO_B_VALUE: { int n = vec.size(); for(int i=0; i void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int ) { // init output and b-zero iterators typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread); oit2.GoToBegin(); vnl_vector B(m_NumberOfGradientDirections); // Two cases here: // 1. Gradients specified in multiple images // 'n' iterators for each of the gradient images // 2. Gradients specified in a single multi-component image // one iterator for all gradient directions if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) { // b-zero reference image iterator ImageRegionConstIterator< ReferenceImageType > it(static_cast< ReferenceImageType * >(this->ProcessObject::GetInput(0)), outputRegionForThread); it.GoToBegin(); // fill vector with gradient iterators typedef ImageRegionConstIterator< GradientImageType > GradientIteratorType; std::vector< GradientIteratorType * > gradientItContainer; for( unsigned int i = 1; i<= m_NumberOfGradientDirections; i++ ) { // adapt index in case we have a duplicated shell int index = i; if(m_DirectionsDuplicated) index = i % (m_NumberOfGradientDirections/2); // init and pushback current input image iterator typename GradientImageType::Pointer gradientImagePointer = NULL; // dynamic_cast would be nice, static because of SGI gradientImagePointer = static_cast< GradientImageType * >( this->ProcessObject::GetInput(index) ); GradientIteratorType *git = new GradientIteratorType( gradientImagePointer, outputRegionForThread ); git->GoToBegin(); gradientItContainer.push_back(git); } // Following loop does the actual reconstruction work in each voxel // (Tuch, Q-Ball Reconstruction [1]) while( !it.IsAtEnd() ) { // b-zero reference value ReferencePixelType b0 = it.Get(); // init ODF OdfPixelType odf(0.0); // threshold on reference value to suppress noisy regions if( (b0 != 0) && (b0 >= m_Threshold) ) { // fill array of diffusion measurements for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { GradientPixelType b = gradientItContainer[i]->Get(); B[i] = static_cast(b); ++(*gradientItContainer[i]); } // pre-normalization according to m_NormalizationMethod B = PreNormalize(B); // actual reconstruction odf = ( (*m_ReconstructionMatrix) * B ).data_block(); // post-normalization according to m_NormalizationMethod odf.Normalize(); } else { // in case we fall below threshold, we just increment to next voxel for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { ++(*gradientItContainer[i]); } } // set and increment output iterators oit.Set( odf ); ++oit; oit2.Set( b0 ); ++oit2; ++it; } // clean up for( unsigned int i = 0; i< gradientItContainer.size(); i++ ) { delete gradientItContainer[i]; } } // The gradients are specified in a single multi-component image else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage ) { // init input iterator typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; typedef typename GradientImagesType::PixelType GradientVectorType; typename GradientImagesType::Pointer gradientImagePointer = NULL; // dynamic_cast would be nice, static because of SGI gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); GradientIteratorType git(gradientImagePointer, outputRegionForThread ); git.GoToBegin(); // set of indicies each for the baseline images and gradient images std::vector baselineind; // contains baseline indicies std::vector gradientind; // contains gradient indicies for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() <= 0.0) { baselineind.push_back(gdcit.Index()); } else { gradientind.push_back(gdcit.Index()); } } // in case we have a duplicated shell, we also duplicate or indices if( m_DirectionsDuplicated ) { int gradIndSize = gradientind.size(); for(int i=0; i::AccumulateType b0 = NumericTraits::Zero; for(unsigned int i = 0; i < baselineind.size(); ++i) { b0 += b[baselineind[i]]; } b0 /= this->m_NumberOfBaselineImages; // init resulting ODF OdfPixelType odf(0.0); // threshold on reference value to suppress noisy regions if( (b0 != 0) && (b0 >= m_Threshold) ) { for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ ) { B[i] = static_cast(b[gradientind[i]]); } // pre-normalization according to m_NormalizationMethod B = PreNormalize(B); // actual reconstruction odf = ( (*m_ReconstructionMatrix) * B ).data_block(); // post-normalization according to m_NormalizationMethod odf = Normalize(odf, b0); } // set and increment output iterators oit.Set( odf ); ++oit; oit2.Set( b0 ); ++oit2; ++git; // Gradient image iterator } } std::cout << "One Thread finished reconstruction" << std::endl; } //void tofile(vnl_matrix A, std::string fname) //{ // ofstream myfile; // myfile.open (fname.c_str()); // myfile << "A1=["; // for(int i=0; i void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::ComputeReconstructionMatrix() { if( m_NumberOfGradientDirections < 6 ) { itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" ); } { // check for duplicate diffusion gradients bool warning = false; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin(); gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2) { if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index()) { itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); warning = true; break; } } if (warning) break; } // handle acquisition schemes where only half of the spherical // shell is sampled by the gradient directions. In this case, // each gradient direction is duplicated in negative direction. vnl_vector centerMass(3); centerMass.fill(0.0); int count = 0; for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) { if(gdcit1.Value().one_norm() > 0.0) { centerMass += gdcit1.Value(); count ++; } } centerMass /= count; if(centerMass.two_norm() > 0.1) { m_DirectionsDuplicated = true; m_NumberOfGradientDirections *= 2; } } // set default number of equator sampling points if needed if(!this->m_NumberOfEquatorSamplingPoints) this->m_NumberOfEquatorSamplingPoints = (int) ceil((double)sqrt(8*QBALL_RECON_PI*this->m_NumberOfGradientDirections)); vnl_matrix* Q = new vnl_matrix(3, m_NumberOfGradientDirections); { // Fill matrix Q with gradient directions int i = 0; for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { (*Q)(0,i) = gdcit.Value().get(0); (*Q)(1,i) = gdcit.Value().get(1); (*Q)(2,i++) = gdcit.Value().get(2); } } if(m_DirectionsDuplicated) { for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) { if(gdcit.Value().one_norm() > 0.0) { (*Q)(0,i) = -gdcit.Value().get(0); (*Q)(1,i) = -gdcit.Value().get(1); (*Q)(2,i++) = -gdcit.Value().get(2); } } } } vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); vnl_matrix_fixed* V = itk::PointShell >::DistributePointShell(); // calculate sampling points on the equator perpendicular to z-axis vnl_matrix *C = new vnl_matrix(3, m_NumberOfEquatorSamplingPoints); for(unsigned int i=0; i::Zero; } // rotate the sampling points to each directions of interest vnl_matrix *S = new vnl_matrix(3,m_NumberOfEquatorSamplingPoints*NOdfDirections); { vnl_vector_fixed z(NumericTraits::Zero); z.put(2,NumericTraits::One); vnl_matrix_fixed eye(NumericTraits::Zero); eye.fill_diagonal(NumericTraits::One); for(int i=0; i ui = (*U).get_column(i); vnl_matrix *RC = new vnl_matrix(3,m_NumberOfEquatorSamplingPoints); if( (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1) != 0 ) { vnl_matrix_fixed R; R.set_column(0, (z+ui)*(z(0)+ui(0))); R.set_column(1, (z+ui)*(z(1)+ui(1))); R.set_column(2, (z+ui)*(z(2)+ui(2))); R /= (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1); R -= eye; (*RC) = R*(*C); } else { RC = C; } (*S).set_columns(i*m_NumberOfEquatorSamplingPoints, *RC); } } // determine interpolation kernel width first // use to determine diffusion measurement contribution to each of the kernels vnl_matrix *H_plus = new vnl_matrix(NBasisFunctionCenters,m_NumberOfGradientDirections); double maxSigma = QBALL_RECON_PI/6; double bestSigma = maxSigma; { double stepsize = 0.01; double start = 0.01; double minCondition = NumericTraits::max(); vnl_matrix *H = new vnl_matrix(m_NumberOfGradientDirections,NBasisFunctionCenters,(double)0); { int increasing = 0; for( double sigma=start; sigma *tmpH = new vnl_matrix(m_NumberOfGradientDirections,NBasisFunctionCenters); for(unsigned int r=0; r1.0) ? 1.0 : qtv); double x = acos(qtv); (*tmpH)(r,c) = (1.0/(sigma*sqrt(2.0*QBALL_RECON_PI))) *exp((-x*x)/(2*sigma*sigma)); } } vnl_svd *solver; if(m_NumberOfGradientDirections>NBasisFunctionCenters) { solver = new vnl_svd(*tmpH); } else { solver = new vnl_svd(tmpH->transpose()); } double condition = solver->sigma_max() / solver->sigma_min(); std::cout << sigma << ": " << condition << std::endl; if( condition < minCondition ) { minCondition = condition; bestSigma = sigma; H->update(*tmpH); } else { // optimum assumed to be hit after condition increased 3 times if (++increasing>3) break; } } } vnl_matrix_inverse *pseudoInverse = new vnl_matrix_inverse(*H); (*H_plus) = pseudoInverse->pinverse(); std::cout << "choosing sigma = " << bestSigma << std::endl; } // this is the contribution of each kernel to each sampling point on the // equator vnl_matrix *G = new vnl_matrix(m_NumberOfEquatorSamplingPoints*NOdfDirections,NBasisFunctionCenters); { for(unsigned int r=0; r1.0) ? 1.0 : stv); double x = acos(stv); (*G)(r,c) = (1.0/(bestSigma*sqrt(2.0*QBALL_RECON_PI))) *exp((-x*x)/(2*bestSigma*bestSigma)); } } } vnl_matrix *GH_plus = new vnl_matrix(m_NumberOfEquatorSamplingPoints*NOdfDirections,m_NumberOfGradientDirections); // simple matrix multiplication, manual cause of stack overflow using operator for (unsigned i = 0; i < m_NumberOfEquatorSamplingPoints*NOdfDirections; ++i) for (unsigned j = 0; j < m_NumberOfGradientDirections; ++j) { double accum = (*G)(i,0) * (*H_plus)(0,j); for (unsigned k = 1; k < NOdfDirections; ++k) accum += (*G)(i,k) * (*H_plus)(k,j); (*GH_plus)(i,j) = accum; } typename vnl_matrix::iterator it3; for( it3 = (*GH_plus).begin(); it3 != (*GH_plus).end(); it3++) { if(*it3<0.0) *it3 = 0; } // this is an addition to the original tuch algorithm for(unsigned int i=0; i r = GH_plus->get_row(i); r /= r.sum(); GH_plus->set_row(i,r); } //tofile((*G),"G.m"); //tofile((*H_plus),"H_plus.m"); //tofile((*GH_plus),"GH_plus.m"); m_ReconstructionMatrix = new vnl_matrix(NOdfDirections,m_NumberOfGradientDirections,0.0); for(int i=0; i r = m_ReconstructionMatrix->get_row(i); r /= r.sum(); m_ReconstructionMatrix->set_row(i,r); } std::cout << "Reconstruction Matrix computed." << std::endl; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::AddGradientImage( const GradientDirectionType &gradientDirection, const GradientImageType *gradientImage ) { // Make sure crazy users did not call both AddGradientImage and // SetGradientImage if( m_GradientImageTypeEnumeration == GradientIsInASingleImage) { itkExceptionMacro( << "Cannot call both methods:" << "AddGradientImage and SetGradientImage. Please call only one of them."); } // If the container to hold the gradient directions hasn't been allocated // yet, allocate it. if( !this->m_GradientDirectionContainer ) { this->m_GradientDirectionContainer = GradientDirectionContainerType::New(); } this->m_NumberOfGradientDirections = m_GradientDirectionContainer->Size(); m_GradientDirectionContainer->InsertElement( this->m_NumberOfGradientDirections, gradientDirection / gradientDirection.two_norm() ); this->ProcessObject::SetNthInput( this->m_NumberOfGradientDirections, const_cast< GradientImageType* >(gradientImage) ); m_GradientImageTypeEnumeration = GradientIsInManyImages; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { // Make sure crazy users did not call both AddGradientImage and // SetGradientImage if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) { itkExceptionMacro( << "Cannot call both methods:" << "AddGradientImage and SetGradientImage. Please call only one of them."); } this->m_GradientDirectionContainer = gradientDirection; unsigned int numImages = gradientDirection->Size(); this->m_NumberOfBaselineImages = 0; for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin(); it != this->m_GradientDirectionContainer->End(); it++) { if(it.Value().one_norm() <= 0.0) { this->m_NumberOfBaselineImages++; } else // Normalize non-zero gradient directions { it.Value() = it.Value() / it.Value().two_norm(); } } this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; // ensure that the gradient image we received has as many components as // the number of gradient directions if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) { itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages << " directions specified but image has " << gradientImage->GetVectorLength() << " components."); } this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); m_GradientImageTypeEnumeration = GradientIsInASingleImage; } template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NrOdfDirections, int NrBasisFunctionCenters> void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters> ::PrintSelf(std::ostream& os, Indent indent) const { + std::locale C("C"); + std::locale originalLocale = os.getloc(); + os.imbue(C); + Superclass::PrintSelf(os,indent); os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl; if ( m_GradientDirectionContainer ) { os << indent << "GradientDirectionContainer: " << m_GradientDirectionContainer << std::endl; } else { os << indent << "GradientDirectionContainer: (Gradient directions not set)" << std::endl; } os << indent << "NumberOfGradientDirections: " << m_NumberOfGradientDirections << std::endl; os << indent << "NumberOfBaselineImages: " << m_NumberOfBaselineImages << std::endl; os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl; os << indent << "BValue: " << m_BValue << std::endl; if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages ) { os << indent << "Gradient images haven been supplied " << std::endl; } else if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages ) { os << indent << "A multicomponent gradient image has been supplied" << std::endl; } + + os.imbue( originalLocale ); } } #endif // __itkDiffusionQballReconstructionImageFilter_txx diff --git a/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp index f1a9827bfb..261d6c0c78 100644 --- a/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp @@ -1,217 +1,228 @@ /*========================================================================= 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 #include #include #include #include #include "mitkPixelType.h" #include "itkImageRegionIterator.h" #include "itkImageFileReader.h" #include "mitkNrrdDiffusionImageWriter.h" template< class D, class T > mitk::TeemDiffusionTensor3DReconstructionImageFilter ::TeemDiffusionTensor3DReconstructionImageFilter(): m_EstimateErrorImage(false),m_Sigma(-19191919), m_EstimationMethod(TeemTensorEstimationMethodsLLS), m_NumIterations(1),m_ConfidenceThreshold(-19191919.0), m_ConfidenceFuzzyness(0.0),m_MinPlausibleValue(1.0) { } template< class D, class T > mitk::TeemDiffusionTensor3DReconstructionImageFilter ::~TeemDiffusionTensor3DReconstructionImageFilter() { } void file_replace(std::string filename, std::string what, std::string with) { ofstream myfile2; + + std::locale C("C"); + std::locale originalLocale2 = myfile2.getloc(); + myfile2.imbue(C); + char filename2[512]; sprintf(filename2, "%s2",filename.c_str()); myfile2.open (filename2); std::string line; ifstream myfile (filename.c_str()); + + std::locale originalLocale = myfile.getloc(); + myfile.imbue(C); + if (myfile.is_open()) { while (! myfile.eof() ) { getline (myfile,line); itksys::SystemTools::ReplaceString(line,what.c_str(),with.c_str()); myfile2 << line << std::endl; } myfile.close(); } myfile2.close(); itksys::SystemTools::RemoveFile(filename.c_str()); rename(filename2,filename.c_str()); + myfile.imbue( originalLocale ); + myfile2.imbue( originalLocale2 ); } // do the work template< class D, class T > void mitk::TeemDiffusionTensor3DReconstructionImageFilter ::Update() { // save input image to nrrd file in temp-folder char filename[512]; srand((unsigned)time(0)); int random_integer = rand(); sprintf( filename, "dwi_%d.nhdr",random_integer); typedef mitk::NrrdDiffusionImageWriter WriterType; typename WriterType::Pointer nrrdWriter = WriterType::New(); nrrdWriter->SetInput( m_Input ); //nrrdWriter->SetDirections(m_Input->GetDirections()); //nrrdWriter->SetB_Value(m_Input->GetB_Value()); nrrdWriter->SetFileName(filename); try { nrrdWriter->Update(); } catch (itk::ExceptionObject e) { std::cout << e << std::endl; } file_replace(filename,"vector","list"); // build up correct command from input params char command[4096]; sprintf( command, "tend estim -i %s -B kvp -o tensors_%d.nhdr -knownB0 true", filename, random_integer); //m_DiffusionImages; if(m_EstimateErrorImage) { sprintf( command, "%s -ee error_image_%d.nhdr", command, random_integer); } if(m_Sigma != -19191919) { sprintf( command, "%s -sigma %f", command, m_Sigma); } switch(m_EstimationMethod) { case TeemTensorEstimationMethodsLLS: sprintf( command, "%s -est lls", command); break; case TeemTensorEstimationMethodsMLE: sprintf( command, "%s -est mle", command); break; case TeemTensorEstimationMethodsNLS: sprintf( command, "%s -est nls", command); break; case TeemTensorEstimationMethodsWLS: sprintf( command, "%s -est wls", command); break; } sprintf( command, "%s -wlsi %d", command, m_NumIterations); if(m_ConfidenceThreshold != -19191919.0) { sprintf( command, "%s -t %f", command, m_ConfidenceThreshold); } sprintf( command, "%s -soft %f", command, m_ConfidenceFuzzyness); sprintf( command, "%s -mv %f", command, m_MinPlausibleValue); // call tend estim command std::cout << "Calling <" << command << ">" << std::endl; int success = system(command); if(!success) { MITK_ERROR << "system command could not be called!"; } remove(filename); sprintf( filename, "dwi_%d.raw", random_integer); remove(filename); // change kind from tensor to vector sprintf( filename, "tensors_%d.nhdr", random_integer); file_replace(filename,"3D-masked-symmetric-matrix","vector"); //file_replace(filename,"3D-symmetric-matrix","vector"); //itksys::SystemTools::ReplaceString(line,"3D-masked-symmetric-matrix","vector"); //itksys::SystemTools::ReplaceString(line,"vector","domain"); // read result as mitk::Image and provide it in m_Output typedef itk::ImageFileReader FileReaderType; typename FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(filename); reader->Update(); typename VectorImageType::Pointer vecImage = reader->GetOutput(); remove(filename); sprintf( filename, "tensors_%d.raw", random_integer); remove(filename); typename ItkTensorImageType::Pointer itkTensorImage = ItkTensorImageType::New(); itkTensorImage->SetSpacing( vecImage->GetSpacing() ); // Set the image spacing itkTensorImage->SetOrigin( vecImage->GetOrigin() ); // Set the image origin itkTensorImage->SetDirection( vecImage->GetDirection() ); // Set the image direction itkTensorImage->SetLargestPossibleRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->SetBufferedRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->SetRequestedRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->Allocate(); itk::ImageRegionIterator it(vecImage, vecImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2(itkTensorImage, itkTensorImage->GetLargestPossibleRegion()); it2 = it2.Begin(); //#pragma omp parallel private (it) { for(it=it.Begin();!it.IsAtEnd(); ++it, ++it2) { //#pragma omp single nowait { VectorType vec = it.Get(); TensorType tensor; for(int i=1;i<7;i++) tensor[i-1] = vec[i] * vec[0]; it2.Set( tensor ); } } // end for } // end ompparallel m_OutputItk = mitk::TensorImage::New(); m_OutputItk->InitializeByItk(itkTensorImage.GetPointer()); m_OutputItk->SetVolume( itkTensorImage->GetBufferPointer() ); // in case: read resulting error-image and provide it in m_ErrorImage if(m_EstimateErrorImage) { // open error image here } }