diff --git a/Examples/Plugins/org.mitk.example.gui.pcaexample/src/internal/PCAExample.cpp b/Examples/Plugins/org.mitk.example.gui.pcaexample/src/internal/PCAExample.cpp index c0c6fa8738..cfeeff9224 100644 --- a/Examples/Plugins/org.mitk.example.gui.pcaexample/src/internal/PCAExample.cpp +++ b/Examples/Plugins/org.mitk.example.gui.pcaexample/src/internal/PCAExample.cpp @@ -1,167 +1,171 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "PCAExample.h" // Qt #include // mitk image #include const std::string PCAExample::VIEW_ID = "org.mitk.views.pcaexample"; void PCAExample::SetFocus() { m_Controls.buttonPerformImageProcessing->setFocus(); } void PCAExample::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect(m_Controls.buttonPerformImageProcessing, SIGNAL(clicked()), this, SLOT(BtnPerfomPCAClicked())); //initialize point set widget and point set node mitk::DataNode::Pointer PointSetNode = mitk::DataNode::New(); PointSetNode->SetName("PCA Example Pointset"); mitk::PointSet::Pointer newPtSet = mitk::PointSet::New(); PointSetNode->SetData(newPtSet); m_Controls.m_pointSetWidget->SetPointSetNode(PointSetNode); this->GetDataStorage()->Add(PointSetNode); } PCAExample::PCAExample() { } PCAExample::~PCAExample() { //clean up mitk::DataNode::Pointer ptSetNode = m_Controls.m_pointSetWidget->GetPointSetNode(); m_Controls.m_pointSetWidget->SetPointSetNode(nullptr); this->GetDataStorage()->Remove(ptSetNode); this->GetDataStorage()->Remove(m_Axis1Node); this->GetDataStorage()->Remove(m_Axis2Node); this->GetDataStorage()->Remove(m_Axis3Node); } void PCAExample::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/, const QList &nodes) { // iterate all selected objects, adjust warning visibility foreach (mitk::DataNode::Pointer node, nodes) { if (node.IsNotNull() && dynamic_cast(node->GetData())) { m_Controls.buttonPerformImageProcessing->setEnabled(true); return; } } } void PCAExample::BtnPerfomPCAClicked() { std::vector eigenVectors; std::vector eigenValues; mitk::Vector3D mean; bool success = comutePCA(m_Controls.m_pointSetWidget->GetPointSet(), eigenVectors, eigenValues, mean); this->showEigenvectors(eigenVectors, eigenValues, mean); MITK_INFO << "PCA: " << success; } bool PCAExample::comutePCA(mitk::PointSet::Pointer input, std::vector &eigenVectors, std::vector &eigenValues, mitk::Vector3D &pointsMean) { //Step 1: Construct data matrix vnl_matrix dataMatrix(3, input->GetSize(), 0.0); int size = input->GetSize(); for (int i=0; iGetPoint(i)[0]; dataMatrix[1][i] = input->GetPoint(i)[1]; dataMatrix[2][i] = input->GetPoint(i)[2]; } //Step 2: Remove average for each row (Mittelwertbefreiung) -mitk::Vector3D mean; +mitk::Vector3D mean; for (int i = 0; i < size; i++) { mean += mitk::Vector3D(dataMatrix.get_column(i)); } mean /= size; -for (int i = 0; i covMatrix = (1.0 / (size - 1.0)) * dataMatrix * dataMatrix.transpose(); //Step 4: Singular value composition vnl_svd svd(covMatrix); //Store results and print them to the console MITK_INFO << "DataMatrix: " << "\n" << dataMatrix; MITK_INFO << "CovMatrix: " << "\n" << covMatrix; for (int i = 0; i < 3; i++) { eigenVectors.push_back(svd.U().get_column(i)); eigenValues.push_back(sqrt(svd.W(i))); MITK_INFO << "Eigenvector " << i << ": " << eigenVectors.at(i); MITK_INFO << "Eigenvalue " << i << ": " << eigenValues.at(i); } //Compute center of points for (int i = 0; i < size; i++) { pointsMean += input->GetPoint(i).GetVectorFromOrigin(); } pointsMean /= size; return true; } void PCAExample::showEigenvectors(std::vector eigenVectors, std::vector eigenValues, mitk::Vector3D center) { m_Axis1Node = mitk::DataNode::New(); m_Axis1Node->SetName("Eigenvector 1"); mitk::PointSet::Pointer axis1 = mitk::PointSet::New(); axis1->InsertPoint(0, center); axis1->InsertPoint(1, (center + eigenVectors.at(0)*eigenValues.at(0))); m_Axis1Node->SetData(axis1); m_Axis1Node->SetBoolProperty("show contour", true); m_Axis1Node->SetColor(1, 0, 0); this->GetDataStorage()->Add(m_Axis1Node); m_Axis2Node = mitk::DataNode::New(); m_Axis2Node->SetName("Eigenvector 2"); mitk::PointSet::Pointer axis2 = mitk::PointSet::New(); axis2->InsertPoint(0, center); axis2->InsertPoint(1, (center + eigenVectors.at(1)*eigenValues.at(1))); m_Axis2Node->SetData(axis2); m_Axis2Node->SetBoolProperty("show contour", true); m_Axis2Node->SetColor(1, 0, 0); this->GetDataStorage()->Add(m_Axis2Node); m_Axis3Node = mitk::DataNode::New(); m_Axis3Node->SetName("Eigenvector 3"); mitk::PointSet::Pointer axis3 = mitk::PointSet::New(); axis3->InsertPoint(0, center); axis3->InsertPoint(1, (center + eigenVectors.at(2)*eigenValues.at(2))); m_Axis3Node->SetData(axis3); m_Axis3Node->SetBoolProperty("show contour", true); m_Axis3Node->SetColor(1, 0, 0); this->GetDataStorage()->Add(m_Axis3Node); } diff --git a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp index 6368711f0a..87d46e5e63 100644 --- a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp @@ -1,343 +1,343 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDiffusionCoreIOMimeTypes.h" #include "mitkIOMimeTypes.h" #include #include #include #include #include namespace mitk { std::vector DiffusionCoreIOMimeTypes::Get() { std::vector mimeTypes; // order matters here (descending rank for mime types) mimeTypes.push_back(DWI_NRRD_MIMETYPE().Clone()); mimeTypes.push_back(DWI_NIFTI_MIMETYPE().Clone()); mimeTypes.push_back(DWI_FSL_MIMETYPE().Clone()); mimeTypes.push_back(DTI_MIMETYPE().Clone()); mimeTypes.push_back(QBI_MIMETYPE().Clone()); return mimeTypes; } // Mime Types DiffusionCoreIOMimeTypes::DiffusionImageNrrdMimeType::DiffusionImageNrrdMimeType() : CustomMimeType(DWI_NRRD_MIMETYPE_NAME()) { std::string category = "Diffusion Weighted Images"; this->SetCategory(category); this->SetComment("Diffusion Weighted Images"); this->AddExtension("dwi"); //this->AddExtension("hdwi"); // saving with detached header does not work out of the box this->AddExtension("nrrd"); } bool DiffusionCoreIOMimeTypes::DiffusionImageNrrdMimeType::AppliesTo(const std::string &path) const { bool canRead( CustomMimeType::AppliesTo(path) ); // fix for bug 18572 // Currently this function is called for writing as well as reading, in that case // the image information can of course not be read // This is a bug, this function should only be called for reading. if( ! itksys::SystemTools::FileExists( path.c_str() ) ) { return canRead; } //end fix for bug 18572 std::string ext = this->GetExtension( path ); ext = itksys::SystemTools::LowerCase( ext ); // Simple NRRD files should only be considered for this mime type if they contain // corresponding tags if( ext == ".nrrd" ) { itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); io->SetFileName(path); try { io->ReadImageInformation(); itk::MetaDataDictionary imgMetaDictionary = io->GetMetaDataDictionary(); std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); std::vector::const_iterator itKey = imgMetaKeys.begin(); std::string metaString; for (; itKey != imgMetaKeys.end(); itKey ++) { itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); if (itKey->find("modality") != std::string::npos) { if (metaString.find("DWMRI") != std::string::npos) { return canRead; } } } } catch( const itk::ExceptionObject &e ) { MITK_ERROR << "ITK Exception: " << e.what(); } canRead = false; } return canRead; } DiffusionCoreIOMimeTypes::DiffusionImageNrrdMimeType* DiffusionCoreIOMimeTypes::DiffusionImageNrrdMimeType::Clone() const { return new DiffusionImageNrrdMimeType(*this); } DiffusionCoreIOMimeTypes::DiffusionImageNrrdMimeType DiffusionCoreIOMimeTypes::DWI_NRRD_MIMETYPE() { return DiffusionImageNrrdMimeType(); } DiffusionCoreIOMimeTypes::DiffusionImageNiftiMimeType::DiffusionImageNiftiMimeType() : CustomMimeType(DWI_NIFTI_MIMETYPE_NAME()) { std::string category = "Diffusion Weighted Images"; this->SetCategory(category); this->SetComment("Diffusion Weighted Images"); this->AddExtension("nii.gz"); this->AddExtension("nii"); } bool DiffusionCoreIOMimeTypes::DiffusionImageNiftiMimeType::AppliesTo(const std::string &path) const { MITK_INFO << "path1: " << path; bool canRead(CustomMimeType::AppliesTo(path)); // fix for bug 18572 // Currently this function is called for writing as well as reading, in that case // the image information can of course not be read // This is a bug, this function should only be called for reading. if (!itksys::SystemTools::FileExists(path.c_str())) { return canRead; } //end fix for bug 18572 std::string ext = this->GetExtension(path); ext = itksys::SystemTools::LowerCase(ext); // Nifti files should only be considered for this mime type if they are // accompanied by bvecs and bvals files defining the diffusion information if (ext == ".nii" || ext == ".nii.gz") { std::string base_path = itksys::SystemTools::GetFilenamePath(path); std::string base = this->GetFilenameWithoutExtension(path); if (!base_path.empty()) - std::string base = base_path + "/" + base; + base = base_path + "/" + base; if (itksys::SystemTools::FileExists(std::string(base + ".bvec").c_str()) && itksys::SystemTools::FileExists(std::string(base + ".bval").c_str()) ) { return canRead; } if (itksys::SystemTools::FileExists(std::string(base + ".bvecs").c_str()) && itksys::SystemTools::FileExists(std::string(base + ".bvals").c_str()) ) { return canRead; } canRead = false; } return canRead; } DiffusionCoreIOMimeTypes::DiffusionImageNiftiMimeType* DiffusionCoreIOMimeTypes::DiffusionImageNiftiMimeType::Clone() const { return new DiffusionImageNiftiMimeType(*this); } DiffusionCoreIOMimeTypes::DiffusionImageNiftiMimeType DiffusionCoreIOMimeTypes::DWI_NIFTI_MIMETYPE() { return DiffusionImageNiftiMimeType(); } DiffusionCoreIOMimeTypes::DiffusionImageFslMimeType::DiffusionImageFslMimeType() : CustomMimeType(DWI_FSL_MIMETYPE_NAME()) { std::string category = "Diffusion Weighted Image"; this->SetCategory(category); this->SetComment("Diffusion Weighted Images"); this->AddExtension("fslgz"); this->AddExtension("fsl"); } bool DiffusionCoreIOMimeTypes::DiffusionImageFslMimeType::AppliesTo(const std::string &path) const { bool canRead(CustomMimeType::AppliesTo(path)); // fix for bug 18572 // Currently this function is called for writing as well as reading, in that case // the image information can of course not be read // This is a bug, this function should only be called for reading. if (!itksys::SystemTools::FileExists(path.c_str())) { return canRead; } //end fix for bug 18572 std::string ext = this->GetExtension(path); ext = itksys::SystemTools::LowerCase(ext); // Nifti files should only be considered for this mime type if they are // accompanied by bvecs and bvals files defining the diffusion information if (ext == ".fsl" || ext == ".fslgz") { std::string base_path = itksys::SystemTools::GetFilenamePath(path); std::string base = this->GetFilenameWithoutExtension(path); if (!base_path.empty()) - std::string base = base_path + "/" + base; + base = base_path + "/" + base; if (itksys::SystemTools::FileExists(std::string(base + ".bvec").c_str()) && itksys::SystemTools::FileExists(std::string(base + ".bval").c_str()) ) { return canRead; } if (itksys::SystemTools::FileExists(std::string(base + ".bvecs").c_str()) && itksys::SystemTools::FileExists(std::string(base + ".bvals").c_str()) ) { return canRead; } canRead = false; } return canRead; } DiffusionCoreIOMimeTypes::DiffusionImageFslMimeType* DiffusionCoreIOMimeTypes::DiffusionImageFslMimeType::Clone() const { return new DiffusionImageFslMimeType(*this); } DiffusionCoreIOMimeTypes::DiffusionImageFslMimeType DiffusionCoreIOMimeTypes::DWI_FSL_MIMETYPE() { return DiffusionImageFslMimeType(); } CustomMimeType DiffusionCoreIOMimeTypes::DTI_MIMETYPE() { CustomMimeType mimeType(DTI_MIMETYPE_NAME()); std::string category = "Tensor Images"; mimeType.SetComment("Diffusion Tensor Images"); mimeType.SetCategory(category); mimeType.AddExtension("dti"); //mimeType.AddExtension("hdti"); // saving with detached header does not work out of the box return mimeType; } CustomMimeType DiffusionCoreIOMimeTypes::QBI_MIMETYPE() { CustomMimeType mimeType(QBI_MIMETYPE_NAME()); std::string category = "Q-Ball Images"; mimeType.SetComment("Diffusion Q-Ball Images"); mimeType.SetCategory(category); mimeType.AddExtension("qbi"); //mimeType.AddExtension("hqbi"); // saving with detached header does not work out of the box return mimeType; } // Names std::string DiffusionCoreIOMimeTypes::DWI_NRRD_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".dwi"; return name; } std::string DiffusionCoreIOMimeTypes::DWI_NIFTI_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".nii.gz"; return name; } std::string DiffusionCoreIOMimeTypes::DWI_FSL_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".fslgz"; return name; } std::string DiffusionCoreIOMimeTypes::DTI_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".dti"; return name; } std::string DiffusionCoreIOMimeTypes::QBI_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".qbi"; return name; } // Descriptions std::string DiffusionCoreIOMimeTypes::DWI_NRRD_MIMETYPE_DESCRIPTION() { static std::string description = "Diffusion Weighted Images"; return description; } std::string DiffusionCoreIOMimeTypes::DWI_NIFTI_MIMETYPE_DESCRIPTION() { static std::string description = "Diffusion Weighted Images"; return description; } std::string DiffusionCoreIOMimeTypes::DWI_FSL_MIMETYPE_DESCRIPTION() { static std::string description = "Diffusion Weighted Images"; return description; } std::string DiffusionCoreIOMimeTypes::DTI_MIMETYPE_DESCRIPTION() { static std::string description = "Diffusion Tensor Images"; return description; } std::string DiffusionCoreIOMimeTypes::QBI_MIMETYPE_DESCRIPTION() { static std::string description = "Q-Ball Images"; return description; } } diff --git a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionImageNiftiReaderService.cpp b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionImageNiftiReaderService.cpp index 53f6df7746..da6e65aad8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionImageNiftiReaderService.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionImageNiftiReaderService.cpp @@ -1,456 +1,456 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkDiffusionImageNiftiReaderService_cpp #define __mitkDiffusionImageNiftiReaderService_cpp #include "mitkDiffusionImageNiftiReaderService.h" #include #include // Diffusion properties #include #include #include #include // ITK includes #include #include #include "itksys/SystemTools.hxx" #include "itkImageFileReader.h" #include "itkMetaDataObject.h" #include "itkNiftiImageIO.h" #include "mitkCustomMimeType.h" #include "mitkDiffusionCoreIOMimeTypes.h" #include #include #include #include #include "mitkIOUtil.h" namespace mitk { DiffusionImageNiftiReaderService:: DiffusionImageNiftiReaderService(const DiffusionImageNiftiReaderService & other) : AbstractFileReader(other) { } DiffusionImageNiftiReaderService* DiffusionImageNiftiReaderService::Clone() const { return new DiffusionImageNiftiReaderService(*this); } DiffusionImageNiftiReaderService:: ~DiffusionImageNiftiReaderService() {} DiffusionImageNiftiReaderService:: DiffusionImageNiftiReaderService(CustomMimeType mime_type, std::string mime_type_description ) : mitk::AbstractFileReader( mime_type, mime_type_description ) { m_ServiceReg = this->RegisterService(); } //DiffusionImageNiftiReaderService:: //DiffusionImageNiftiReaderService() : mitk::AbstractFileReader( CustomMimeType( mitk::DiffusionCoreIOMimeTypes::DWI_NIFTI_MIMETYPE() ), mitk::DiffusionCoreIOMimeTypes::DWI_NIFTI_MIMETYPE_DESCRIPTION() ) //{ // m_ServiceReg = this->RegisterService(); //} std::vector > DiffusionImageNiftiReaderService:: Read() { std::vector > result; // Since everything is completely read in GenerateOutputInformation() it is stored // in a cache variable. A timestamp is associated. // If the timestamp of the cache variable is newer than the MTime, we only need to // assign the cache variable to the DataObject. // Otherwise, the tree must be read again from the file and OuputInformation must // be updated! if(m_OutputCache.IsNull()) InternalRead(); result.push_back(m_OutputCache.GetPointer()); return result; } void DiffusionImageNiftiReaderService::InternalRead() { OutputType::Pointer outputForCache = OutputType::New(); if ( this->GetInputLocation() == "") { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename to be read is empty!"); } else { try { mitk::LocaleSwitch localeSwitch("C"); MITK_INFO << "DiffusionImageNiftiReaderService: reading image information"; VectorImageType::Pointer itkVectorImage; std::string ext = this->GetMimeType()->GetExtension( this->GetInputLocation() ); ext = itksys::SystemTools::LowerCase( ext ); if(ext == ".fsl" || ext == ".fslgz") { // create temporary file with correct ending for nifti-io std::string fname3 = "temp_dwi"; fname3 += ext == ".fsl" ? ".nii" : ".nii.gz"; itksys::SystemTools::CopyAFile(this->GetInputLocation().c_str(), fname3.c_str()); // create reader and read file typedef itk::Image ImageType4D; itk::NiftiImageIO::Pointer io2 = itk::NiftiImageIO::New(); typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(fname3); reader->SetImageIO(io2); reader->Update(); ImageType4D::Pointer img4 = reader->GetOutput(); // delete temporary file itksys::SystemTools::RemoveFile(fname3.c_str()); // convert 4D file to vector image itkVectorImage = VectorImageType::New(); VectorImageType::SpacingType spacing; ImageType4D::SpacingType spacing4 = img4->GetSpacing(); for(int i=0; i<3; i++) spacing[i] = spacing4[i]; itkVectorImage->SetSpacing( spacing ); // Set the image spacing VectorImageType::PointType origin; ImageType4D::PointType origin4 = img4->GetOrigin(); for(int i=0; i<3; i++) origin[i] = origin4[i]; itkVectorImage->SetOrigin( origin ); // Set the image origin VectorImageType::DirectionType direction; ImageType4D::DirectionType direction4 = img4->GetDirection(); for(int i=0; i<3; i++) for(int j=0; j<3; j++) direction[i][j] = direction4[i][j]; itkVectorImage->SetDirection( direction ); // Set the image direction VectorImageType::RegionType region; ImageType4D::RegionType region4 = img4->GetLargestPossibleRegion(); VectorImageType::RegionType::SizeType size; ImageType4D::RegionType::SizeType size4 = region4.GetSize(); for(int i=0; i<3; i++) size[i] = size4[i]; VectorImageType::RegionType::IndexType index; ImageType4D::RegionType::IndexType index4 = region4.GetIndex(); for(int i=0; i<3; i++) index[i] = index4[i]; region.SetSize(size); region.SetIndex(index); itkVectorImage->SetRegions( region ); itkVectorImage->SetVectorLength(size4[3]); itkVectorImage->Allocate(); itk::ImageRegionIterator it ( itkVectorImage, itkVectorImage->GetLargestPossibleRegion() ); typedef VectorImageType::PixelType VecPixType; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { VecPixType vec = it.Get(); VectorImageType::IndexType currentIndex = it.GetIndex(); for(int i=0; i<3; i++) index4[i] = currentIndex[i]; for(unsigned int ind=0; indGetPixel(index4); } it.Set(vec); } } else if(ext == ".nii" || ext == ".nii.gz") { // create reader and read file typedef itk::Image ImageType4D; itk::NiftiImageIO::Pointer io2 = itk::NiftiImageIO::New(); typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName( this->GetInputLocation() ); reader->SetImageIO(io2); reader->Update(); ImageType4D::Pointer img4 = reader->GetOutput(); // convert 4D file to vector image itkVectorImage = VectorImageType::New(); VectorImageType::SpacingType spacing; ImageType4D::SpacingType spacing4 = img4->GetSpacing(); for(int i=0; i<3; i++) spacing[i] = spacing4[i]; itkVectorImage->SetSpacing( spacing ); // Set the image spacing VectorImageType::PointType origin; ImageType4D::PointType origin4 = img4->GetOrigin(); for(int i=0; i<3; i++) origin[i] = origin4[i]; itkVectorImage->SetOrigin( origin ); // Set the image origin VectorImageType::DirectionType direction; ImageType4D::DirectionType direction4 = img4->GetDirection(); for(int i=0; i<3; i++) for(int j=0; j<3; j++) direction[i][j] = direction4[i][j]; itkVectorImage->SetDirection( direction ); // Set the image direction VectorImageType::RegionType region; ImageType4D::RegionType region4 = img4->GetLargestPossibleRegion(); VectorImageType::RegionType::SizeType size; ImageType4D::RegionType::SizeType size4 = region4.GetSize(); for(int i=0; i<3; i++) size[i] = size4[i]; VectorImageType::RegionType::IndexType index; ImageType4D::RegionType::IndexType index4 = region4.GetIndex(); for(int i=0; i<3; i++) index[i] = index4[i]; region.SetSize(size); region.SetIndex(index); itkVectorImage->SetRegions( region ); itkVectorImage->SetVectorLength(size4[3]); itkVectorImage->Allocate(); itk::ImageRegionIterator it ( itkVectorImage, itkVectorImage->GetLargestPossibleRegion() ); typedef VectorImageType::PixelType VecPixType; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { VecPixType vec = it.Get(); VectorImageType::IndexType currentIndex = it.GetIndex(); for(int i=0; i<3; i++) index4[i] = currentIndex[i]; for(unsigned int ind=0; indGetPixel(index4); } it.Set(vec); } } // Diffusion Image information START GradientDirectionContainerType::Pointer DiffusionVectors = GradientDirectionContainerType::New(); MeasurementFrameType MeasurementFrame; float BValue = -1; // Diffusion Image information END if(ext == ".fsl" || ext == ".fslgz" || ext == ".nii" || ext == ".nii.gz") { // some parsing depending on the extension bool useFSLstyle( true ); std::string bvecsExtension(""); std::string bvalsExtension(""); std::string base_path = itksys::SystemTools::GetFilenamePath(this->GetInputLocation()); std::string base = this->GetMimeType()->GetFilenameWithoutExtension(this->GetInputLocation()); if (!base_path.empty()) - std::string base = base_path + "/" + base; + base = base_path + "/" + base; // check for possible file names { if( useFSLstyle && itksys::SystemTools::FileExists( std::string( base + ".bvec").c_str() ) && itksys::SystemTools::FileExists( std::string( base + ".bval").c_str() ) ) { useFSLstyle = false; bvecsExtension = ".bvec"; bvalsExtension = ".bval"; } if( useFSLstyle && itksys::SystemTools::FileExists( std::string( base + ".bvecs").c_str() ) && itksys::SystemTools::FileExists( std::string( base + ".bvals").c_str() ) ) { useFSLstyle = false; bvecsExtension = ".bvecs"; bvalsExtension = ".bvals"; } } std::string line; std::vector bvec_entries; std::string fname = this->GetInputLocation(); if( useFSLstyle ) { fname += ".bvecs"; } else { fname = std::string( base + bvecsExtension); } std::ifstream myfile (fname.c_str()); if (myfile.is_open()) { while ( myfile.good() ) { std::string line2; getline (myfile,line2); std::stringstream iss; iss << line2; while(getline(iss,line, ' ')) { // remove any potenial control sequences that might be introduced by lines ending in a single space line.erase(std::remove_if(line.begin(), line.end(), [](char c) { return std::isspace(c) || std::iscntrl(c); } ), line.end()); if (line.length() > 0 ) // otherwise string contained only control characters before, empty string are converted to 0 by atof resulting in a broken b-value list { bvec_entries.push_back(atof(line.c_str())); } } } myfile.close(); } else { MITK_INFO << "Unable to open bvecs file. Expected name: " << fname; } std::vector bval_entries; std::string fname2 = this->GetInputLocation(); if( useFSLstyle ) { fname2 += ".bvals"; } else { fname2 = std::string( base + bvalsExtension); } std::ifstream myfile2 (fname2.c_str()); if (myfile2.is_open()) { while ( myfile2.good() ) { getline (myfile2,line, ' '); // remove any potenial control sequences that might be introduced by lines ending in a single space line.erase(std::remove_if(line.begin(), line.end(), [](char c) { return std::isspace(c) || std::iscntrl(c); } ), line.end()); if (line.length() > 0 ) // otherwise string contained only control characters before, empty string are converted to 0 by atof resulting in a broken b-value list { bval_entries.push_back(atof(line.c_str())); } } myfile2.close(); } else { MITK_INFO << "Unable to open bvals file. Expected name: " << fname2; } BValue = -1; unsigned int numb = bval_entries.size(); for(unsigned int i=0; i vec; vec[0] = bvec_entries.at(i); vec[1] = bvec_entries.at(i+numb); vec[2] = bvec_entries.at(i+2*numb); // Adjust the vector length to encode gradient strength float factor = b_val/BValue; if(vec.magnitude() > 0) { vec[0] = sqrt(factor)*vec[0]; vec[1] = sqrt(factor)*vec[1]; vec[2] = sqrt(factor)*vec[2]; } DiffusionVectors->InsertElement(i,vec); } for(int i=0; i<3; i++) for(int j=0; j<3; j++) MeasurementFrame[i][j] = i==j ? 1 : 0; } outputForCache = mitk::GrabItkImageMemory( itkVectorImage); // create BValueMap mitk::BValueMapProperty::BValueMap BValueMap = mitk::BValueMapProperty::CreateBValueMap(DiffusionVectors,BValue); outputForCache->SetProperty( mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( DiffusionVectors ) ); outputForCache->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( MeasurementFrame ) ); outputForCache->SetProperty( mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str(), mitk::BValueMapProperty::New( BValueMap ) ); outputForCache->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( BValue ) ); mitk::DiffusionPropertyHelper propertyHelper( outputForCache ); propertyHelper.InitializeImage(); // Since we have already read the tree, we can store it in a cache variable // so that it can be assigned to the DataObject in GenerateData(); m_OutputCache = outputForCache; m_CacheTime.Modified(); } catch(std::exception& e) { MITK_INFO << "Std::Exception while reading file!!"; MITK_INFO << e.what(); throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what()); } catch(...) { MITK_INFO << "Exception while reading file!!"; throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested vessel tree file!"); } } } } //namespace MITK #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx index 7253e35b42..01798813da 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx @@ -1,261 +1,261 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt 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 _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_txx_ #define _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_txx_ #endif #define _USE_MATH_DEFINES #include "itkElectrostaticRepulsionDiffusionGradientReductionFilter.h" #include #include #include #include namespace itk { template ElectrostaticRepulsionDiffusionGradientReductionFilter ::ElectrostaticRepulsionDiffusionGradientReductionFilter() { this->SetNumberOfRequiredInputs( 1 ); } template double ElectrostaticRepulsionDiffusionGradientReductionFilter ::Costs() { double costs = 2*M_PI; for (IndicesVector::iterator it = m_UsedGradientIndices.begin(); it!=m_UsedGradientIndices.end(); ++it) { for (IndicesVector::iterator it2 = m_UsedGradientIndices.begin(); it2!=m_UsedGradientIndices.end(); ++it2) if (it != it2) { vnl_vector_fixed v1 = m_OriginalGradientDirections->at(*it); vnl_vector_fixed v2 = m_OriginalGradientDirections->at(*it2); v1.normalize(); v2.normalize(); double temp = dot_product(v1,v2); if (temp>1) temp = 1; else if (temp<-1) temp = -1; double angle = acos(temp); if (angle1) temp = 1; else if (temp<-1) temp = -1; angle = acos(temp); if (angle void ElectrostaticRepulsionDiffusionGradientReductionFilter ::GenerateData() { unsigned int randSeed = time(NULL); if(m_InputBValueMap.empty() || m_NumGradientDirections.size()!=m_InputBValueMap.size()) - return; + mitkThrow() << "Vector of the number of desired gradient directions contains more elements than the specified target b-value map."; BValueMap manipulatedMap = m_InputBValueMap; int shellCounter = 0; for(BValueMap::iterator it = m_InputBValueMap.begin(); it != m_InputBValueMap.end(); it++ ) { srand(randSeed); // initialize index vectors m_UsedGradientIndices.clear(); m_UnusedGradientIndices.clear(); if ( it->second.size() <= m_NumGradientDirections[shellCounter] ) { itkWarningMacro( << "current directions: " << it->second.size() << " wanted directions: " << m_NumGradientDirections[shellCounter]); m_NumGradientDirections[shellCounter] = it->second.size(); shellCounter++; continue; } MITK_INFO << "Shell number: " << shellCounter; unsigned int c=0; for (unsigned int i=0; isecond.size(); i++) { if (csecond.at(i)); else m_UnusedGradientIndices.push_back(it->second.at(i)); c++; } double minAngle = Costs(); double newMinAngle = 0; MITK_INFO << "minimum angle: " << 180*minAngle/M_PI; int stagnationCount = 0; int rejectionCount = 0; int maxRejections = m_NumGradientDirections[shellCounter] * 1000; if (maxRejections<10000) maxRejections = 10000; int iUsed = 0; if (m_UsedGradientIndices.size()>0) while ( stagnationCount<1000 && rejectionCount minAngle) // accept or reject proposal { MITK_INFO << "minimum angle: " << 180*newMinAngle/M_PI; if ( (newMinAngle-minAngle)<0.01 ) stagnationCount++; else stagnationCount = 0; minAngle = newMinAngle; rejectionCount = 0; } else { rejectionCount++; m_UsedGradientIndices.at(iUsed) = vUsed; m_UnusedGradientIndices.at(iUnUsed) = vUnUsed; } iUsed++; iUsed = iUsed % m_UsedGradientIndices.size(); } manipulatedMap[it->first] = m_UsedGradientIndices; shellCounter++; } int vecLength = 0 ; for(BValueMap::iterator it = manipulatedMap.begin(); it != manipulatedMap.end(); it++) vecLength += it->second.size(); // initialize output image typename OutputImageType::Pointer outImage = OutputImageType::New(); outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion()); outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetVectorLength( vecLength ); // Set the vector length outImage->Allocate(); itk::ImageRegionIterator< OutputImageType > newIt(outImage, outImage->GetLargestPossibleRegion()); newIt.GoToBegin(); typename InputImageType::Pointer inImage = const_cast(this->GetInput(0)); itk::ImageRegionIterator< InputImageType > oldIt(inImage, inImage->GetLargestPossibleRegion()); oldIt.GoToBegin(); // initial new value of voxel OutputPixelType newVec; newVec.SetSize( vecLength ); newVec.AllocateElements( vecLength ); // generate new pixel values while(!newIt.IsAtEnd()) { // init new vector with zeros newVec.Fill(0.0); InputPixelType oldVec = oldIt.Get(); int index = 0; for(BValueMap::iterator it=manipulatedMap.begin(); it!=manipulatedMap.end(); it++) for(unsigned int j=0; jsecond.size(); j++) { newVec[index] = oldVec[it->second.at(j)]; index++; } newIt.Set(newVec); ++newIt; ++oldIt; } // set new gradient directions m_GradientDirections = GradientDirectionContainerType::New(); int index = 0; for(BValueMap::iterator it = manipulatedMap.begin(); it != manipulatedMap.end(); it++) for(unsigned int j = 0; j < it->second.size(); j++) { m_GradientDirections->InsertElement(index, m_OriginalGradientDirections->at(it->second.at(j))); index++; } this->SetNumberOfRequiredOutputs (1); this->SetNthOutput (0, outImage); MITK_INFO << "...done"; } template void ElectrostaticRepulsionDiffusionGradientReductionFilter ::UpdateOutputInformation() { Superclass::UpdateOutputInformation(); int vecLength = 0 ; for(unsigned int index = 0; index < m_NumGradientDirections.size(); index++) { vecLength += m_NumGradientDirections[index]; } this->GetOutput()->SetVectorLength( vecLength ); } } // end of namespace diff --git a/Modules/DiffusionImaging/cmdapps/CMakeLists.txt b/Modules/DiffusionImaging/cmdapps/CMakeLists.txt index 6ee4d6b149..8e2417bd78 100755 --- a/Modules/DiffusionImaging/cmdapps/CMakeLists.txt +++ b/Modules/DiffusionImaging/cmdapps/CMakeLists.txt @@ -1,59 +1,60 @@ if(BUILD_DiffusionCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusioncmdapps NetworkCreation^^MitkFiberTracking_MitkConnectomics NetworkStatistics^^MitkConnectomics Fiberfox^^MitkFiberTracking MultishellMethods^^MitkFiberTracking PeaksAngularError^^MitkFiberTracking PeakExtraction^^MitkFiberTracking FiberExtraction^^MitkFiberTracking FiberProcessing^^MitkFiberTracking FiberDirectionExtraction^^MitkFiberTracking # LocalDirectionalFiberPlausibility^^MitkFiberTracking # HAS TO USE NEW PEAK IMAGE FORMAT StreamlineTracking^^MitkFiberTracking GibbsTracking^^MitkFiberTracking TractometerMetrics^^MitkFiberTracking FileFormatConverter^^MitkFiberTracking RfTraining^^MitkFiberTracking TractDensity^^MitkFiberTracking Sift2WeightCopy^^MitkFiberTracking + ResampleGradients^^MitkFiberTracking ) foreach(diffusioncmdapp ${diffusioncmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusioncmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} PACKAGE_DEPENDS ITK ) endforeach() # This cmd app does not depend on mitkDiffusionImaging at all mitkFunctionCreateCommandLineApp( NAME Dicom2Nrrd DEPENDS MitkCore ${dependencies_list} ) if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/cmdapps/ResampleGradients.cpp b/Modules/DiffusionImaging/cmdapps/ResampleGradients.cpp new file mode 100644 index 0000000000..d3f7795e70 --- /dev/null +++ b/Modules/DiffusionImaging/cmdapps/ResampleGradients.cpp @@ -0,0 +1,255 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include "mitkCommandLineParser.h" +#include +#include +#include +#include +#include + +#include +#include +#include + + +// itk includes +#include + +// mitk includes +#include +#include "itkDWIVoxelFunctor.h" +#include + +typedef short DiffusionPixelType; +typedef itk::VectorImage< short, 3 > ItkDwiType; + +// itk includes +#include "itkTimeProbe.h" +#include "itkB0ImageExtractionImageFilter.h" +#include "itkB0ImageExtractionToSeparateImageFilter.h" +#include "itkBrainMaskExtractionImageFilter.h" +#include "itkCastImageFilter.h" +#include "itkVectorContainer.h" +#include +#include +#include +#include +#include +#include + +// Multishell includes +#include + +// Multishell Functors +#include +#include +#include +#include + +// mitk includes +#include "mitkProgressBar.h" +#include "mitkStatusBar.h" +#include "mitkNodePredicateDataType.h" +#include "mitkProperties.h" +#include "mitkVtkResliceInterpolationProperty.h" +#include "mitkLookupTable.h" +#include "mitkLookupTableProperty.h" +#include "mitkTransferFunction.h" +#include "mitkTransferFunctionProperty.h" +//#include "mitkDataNodeObject.h" +#include "mitkOdfNormalizationMethodProperty.h" +#include "mitkOdfScaleByProperty.h" +#include +#include +#include +#include +#include +//#include +//#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "mitkPreferenceListReaderOptionsFunctor.h" + +mitk::Image::Pointer DoReduceGradientDirections(mitk::Image::Pointer image, double BValue, unsigned int numOfGradientsToKeep) +{ + + bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(image) ); + if ( !isDiffusionImage ) { + std::cout << "Image is not a Diffusion Weighted Image" << std::endl; + //return; + } + + typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter FilterType; + typedef mitk::BValueMapProperty::BValueMap BValueMap; + + BValueMap shellSlectionMap; + BValueMap originalShellMap = static_cast + (image->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() ) + ->GetBValueMap(); + + std::vector newNumGradientDirections; + + //Keeps 1 b0 gradient + double B0Value = 0; + shellSlectionMap[B0Value] = originalShellMap[B0Value]; + unsigned int num = 1; + newNumGradientDirections.push_back(num); + + //BValue = 1000; + shellSlectionMap[BValue] = originalShellMap[BValue]; + //numOfGradientsToKeep = 32; + newNumGradientDirections.push_back(numOfGradientsToKeep); + + + if (newNumGradientDirections.empty()) { + std::cout << "newNumGradientDirections is empty" << std::endl; + //return; + } + + itk::DwiGradientLengthCorrectionFilter::GradientDirectionContainerType::Pointer gradientContainer + = static_cast + ( image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ) + ->GetGradientDirectionsContainer(); + + + ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); + mitk::CastToItkImage(image, itkVectorImagePointer); + std::cout << "1" << std::endl; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput( itkVectorImagePointer ); + filter->SetOriginalGradientDirections(gradientContainer); + filter->SetNumGradientDirections(newNumGradientDirections); + filter->SetOriginalBValueMap(originalShellMap); + filter->SetShellSelectionBValueMap(shellSlectionMap); + filter->Update(); + std::cout << "2" << std::endl; + + if( filter->GetOutput() == nullptr) { + std::cout << "filter get output is nullptr" << std::endl; + } + + mitk::Image::Pointer newImage = mitk::GrabItkImageMemory( filter->GetOutput() ); + + newImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), + mitk::GradientDirectionsProperty::New( filter->GetGradientDirections() ) ); + + newImage->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), + mitk::MeasurementFrameProperty::New( static_cast + (image->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() ) + ->GetMeasurementFrame() ) ); + + newImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), + mitk::FloatProperty::New( static_cast + (image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() ) + ->GetValue() ) ); + + mitk::DiffusionPropertyHelper propertyHelper( newImage ); + propertyHelper.InitializeImage(); //needed? + + return newImage; +} + + + +/*! +\brief Resample gradients of input DWI image. +*/ +int main(int argc, char* argv[]) +{ + mitkCommandLineParser parser; + + parser.setTitle("Resample Gradients"); + parser.setCategory("Fiber Tracking and Processing Methods"); + parser.setDescription("Resample gradients of input DWI image. You can select one b-value shell and the number of gradients within this shell you want to have. It will also keep one b0 image."); + parser.setContributor("MBI"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle", us::Any(), false); + parser.addArgument("output", "o", mitkCommandLineParser::String, "Output:", "output fiber bundle", us::Any(), false); + parser.addArgument("bValue", "b", mitkCommandLineParser::Float, "b-value:", "float", 1000, false); + parser.addArgument("nrOfGradients", "n", mitkCommandLineParser::Int, "Nr of gradients:", "integer", 32, false); + + + std::map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + + std::string inFileName = us::any_cast(parsedArgs["input"]); + std::string outFileName = us::any_cast(parsedArgs["output"]); + double bValue = us::any_cast(parsedArgs["bValue"]); + unsigned int nrOfGradients = us::any_cast(parsedArgs["nrOfGradients"]); + + try + { + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({ "Diffusion Weighted Images" }, {}); + mitk::Image::Pointer mitkImage = dynamic_cast(mitk::IOUtil::LoadImage(inFileName, &functor).GetPointer()); + mitk::Image::Pointer newImage = DoReduceGradientDirections(mitkImage, bValue, nrOfGradients); + //mitk::IOUtil::SaveImage(newImage, outFileName); //save as dwi image + mitk::IOUtil::Save(newImage, "application/vnd.mitk.nii.gz", outFileName); //save as nifti image + + } + catch (itk::ExceptionObject e) + { + std::cout << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + std::cout << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + std::cout << "ERROR!?!"; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/Modules/DiffusionImaging/cmdapps/Sift2WeightCopy.cpp b/Modules/DiffusionImaging/cmdapps/Sift2WeightCopy.cpp index 4b31ae6801..caeda36330 100644 --- a/Modules/DiffusionImaging/cmdapps/Sift2WeightCopy.cpp +++ b/Modules/DiffusionImaging/cmdapps/Sift2WeightCopy.cpp @@ -1,120 +1,108 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Import Sift2 Fiber Weights txt file. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Sift2 Fiber Weight Import"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Import sift2 fiber weights."); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle", us::Any(), false); parser.addArgument("weights", "w", mitkCommandLineParser::String, "Weights:", "input weights file (.txt)", us::Any(), false); parser.addArgument("output", "o", mitkCommandLineParser::String, "Output:", "output fiber bundle", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; - bool binary = false; - if (parsedArgs.count("binary")) - binary = us::any_cast(parsedArgs["binary"]); - - bool endpoints = false; - if (parsedArgs.count("endpoints")) - endpoints = us::any_cast(parsedArgs["endpoints"]); - - std::string reference_image = ""; - if (parsedArgs.count("reference_image")) - reference_image = us::any_cast(parsedArgs["reference_image"]); - std::string inFileName = us::any_cast(parsedArgs["input"]); std::string weightsFileName = us::any_cast(parsedArgs["weights"]); std::string outFileName = us::any_cast(parsedArgs["output"]); try { mitk::FiberBundle::Pointer fib = LoadFib(inFileName); std::ifstream fin; fin.open(weightsFileName); if (!fin.good()) return 1; // exit if file not found std::vector weights; for (float d; fin >> d; ) { weights.push_back(d); } for(int i = 0; i != weights.size(); i++) { fib->SetFiberWeight(i, weights[i]); } mitk::IOUtil::SaveBaseData(fib.GetPointer(), outFileName ); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; }