diff --git a/CMakeExternals/MITKData.cmake b/CMakeExternals/MITKData.cmake index 601ccdbf9a..46bf725e21 100644 --- a/CMakeExternals/MITKData.cmake +++ b/CMakeExternals/MITKData.cmake @@ -1,39 +1,39 @@ #----------------------------------------------------------------------------- # MITK Data #----------------------------------------------------------------------------- # Sanity checks if(DEFINED MITK_DATA_DIR AND NOT EXISTS ${MITK_DATA_DIR}) message(FATAL_ERROR "MITK_DATA_DIR variable is defined but corresponds to non-existing directory") endif() set(proj MITK-Data) set(proj_DEPENDENCIES) set(MITK-Data_DEPENDS ${proj}) if(BUILD_TESTING) - set(revision_tag fee70ea2) # first 8 characters of hash-tag + set(revision_tag 7968c5c0) # first 8 characters of hash-tag # ^^^^^^^^ these are just to check correct length of hash part ExternalProject_Add(${proj} SOURCE_DIR ${proj} # GIT_REPOSITORY https://phabricator.mitk.org/diffusion/MD/mitk-data.git # GIT_TAG ${revision_tag} URL ${MITK_THIRDPARTY_DOWNLOAD_PREFIX_URL}/mitk-data_${revision_tag}.tar.gz UPDATE_COMMAND "" CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" DEPENDS ${proj_DEPENDENCIES} ) set(MITK_DATA_DIR ${CMAKE_CURRENT_BINARY_DIR}/${proj}) else() mitkMacroEmptyExternalProject(${proj} "${proj_DEPENDENCIES}") endif(BUILD_TESTING) diff --git a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp index 5763b03118..37bb1546d9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp @@ -1,630 +1,642 @@ /*=================================================================== 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 #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(DWI_DICOM_MIMETYPE().Clone()); mimeTypes.push_back(DTI_MIMETYPE().Clone()); mimeTypes.push_back(ODF_MIMETYPE().Clone()); mimeTypes.push_back(PEAK_MIMETYPE().Clone()); mimeTypes.push_back(SH_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 { 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); std::string filename = base; if (!base_path.empty()) { base = base_path + "/" + base; base_path += "/"; } 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; } // hack for HCP data if ( filename=="data" && itksys::SystemTools::FileExists(std::string(base_path + "bvec").c_str()) && itksys::SystemTools::FileExists(std::string(base_path + "bval").c_str()) ) { return canRead; } if ( filename=="data" && itksys::SystemTools::FileExists(std::string(base_path + "bvecs").c_str()) && itksys::SystemTools::FileExists(std::string(base_path + "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 Images"; 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()) 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; } if (itksys::SystemTools::FileExists(std::string(base + ext + ".bvec").c_str()) && itksys::SystemTools::FileExists(std::string(base + ext + ".bval").c_str()) ) { return canRead; } if (itksys::SystemTools::FileExists(std::string(base + ext + ".bvecs").c_str()) && itksys::SystemTools::FileExists(std::string(base + ext + ".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(); } DiffusionCoreIOMimeTypes::DiffusionImageDicomMimeType::DiffusionImageDicomMimeType() : CustomMimeType(DWI_DICOM_MIMETYPE_NAME()) { std::string category = "Diffusion Weighted Images"; this->SetCategory(category); this->SetComment("Diffusion Weighted Images"); this->AddExtension("gdcm"); this->AddExtension("dcm"); this->AddExtension("DCM"); this->AddExtension("dc3"); this->AddExtension("DC3"); this->AddExtension("ima"); this->AddExtension("img"); } bool DiffusionCoreIOMimeTypes::DiffusionImageDicomMimeType::AppliesTo(const std::string &path) const { itk::GDCMImageIO::Pointer gdcmIO = itk::GDCMImageIO::New(); bool canRead = gdcmIO->CanReadFile(path.c_str()); if (!canRead) return canRead; mitk::DICOMDCMTKTagScanner::Pointer scanner = mitk::DICOMDCMTKTagScanner::New(); mitk::DICOMTag ImageTypeTag(0x0008, 0x0008); mitk::DICOMTag SeriesDescriptionTag(0x0008, 0x103E); mitk::StringList relevantFiles; relevantFiles.push_back(path); scanner->AddTag(ImageTypeTag); scanner->AddTag(SeriesDescriptionTag); scanner->SetInputFiles(relevantFiles); scanner->Scan(); mitk::DICOMTagCache::Pointer tagCache = scanner->GetScanCache(); mitk::DICOMImageFrameList imageFrameList = mitk::ConvertToDICOMImageFrameList(tagCache->GetFrameInfoList()); mitk::DICOMImageFrameInfo *firstFrame = imageFrameList.begin()->GetPointer(); std::string byteString = tagCache->GetTagValue(firstFrame, ImageTypeTag).value; if (byteString.empty()) return false; std::string byteString2 = tagCache->GetTagValue(firstFrame, SeriesDescriptionTag).value; if (byteString2.empty()) return false; if (byteString.find("DIFFUSION")==std::string::npos && byteString2.find("diff")==std::string::npos) return false; // if (byteString.find("NONE")==std::string::npos) // return false; return canRead; } DiffusionCoreIOMimeTypes::DiffusionImageDicomMimeType* DiffusionCoreIOMimeTypes::DiffusionImageDicomMimeType::Clone() const { return new DiffusionImageDicomMimeType(*this); } DiffusionCoreIOMimeTypes::DiffusionImageDicomMimeType DiffusionCoreIOMimeTypes::DWI_DICOM_MIMETYPE() { return DiffusionImageDicomMimeType(); } DiffusionCoreIOMimeTypes::PeakImageMimeType::PeakImageMimeType() : CustomMimeType(PEAK_MIMETYPE_NAME()) { std::string category = "Peak Image"; this->SetCategory(category); this->SetComment("Peak Image"); this->AddExtension("nrrd"); this->AddExtension("nii"); this->AddExtension("nii.gz"); + this->AddExtension("peak"); } bool DiffusionCoreIOMimeTypes::PeakImageMimeType::AppliesTo(const std::string &path) const { + std::string ext = itksys::SystemTools::GetFilenameExtension(path); + if (ext==".peak") + return true; + + try { itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); if ( io->CanReadFile( path.c_str() ) ) { io->SetFileName( path.c_str() ); io->ReadImageInformation(); if ( io->GetPixelType() == itk::ImageIOBase::SCALAR && io->GetNumberOfDimensions()==4 && io->GetDimensions(3)%3==0) return true; } } + catch(...) + {} + + try { itk::NiftiImageIO::Pointer io = itk::NiftiImageIO::New(); if ( io->CanReadFile( path.c_str() ) ) { io->SetFileName( path.c_str() ); io->ReadImageInformation(); if ( io->GetPixelType() == itk::ImageIOBase::SCALAR && io->GetNumberOfDimensions()==4 && io->GetDimensions(3)%3==0) return true; } } + catch(...) + {} return false; } DiffusionCoreIOMimeTypes::PeakImageMimeType* DiffusionCoreIOMimeTypes::PeakImageMimeType::Clone() const { return new PeakImageMimeType(*this); } DiffusionCoreIOMimeTypes::PeakImageMimeType DiffusionCoreIOMimeTypes::PEAK_MIMETYPE() { return PeakImageMimeType(); } DiffusionCoreIOMimeTypes::SHImageMimeType::SHImageMimeType() : CustomMimeType(SH_MIMETYPE_NAME()) { std::string category = "SH Image"; this->SetCategory(category); this->SetComment("SH Image"); this->AddExtension("nii.gz"); this->AddExtension("nii"); this->AddExtension("nrrd"); this->AddExtension("shi"); } bool DiffusionCoreIOMimeTypes::SHImageMimeType::AppliesTo(const std::string &path) const { { itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); try { io->SetFileName( path.c_str() ); io->ReadImageInformation(); if ( io->GetPixelType() == itk::ImageIOBase::SCALAR && io->GetNumberOfDimensions()==4) { switch (io->GetDimensions(3)) { case 6: return true; break; case 15: return true; break; case 28: return true; break; case 45: return true; break; case 66: return true; break; case 91: return true; break; default : return false; } } } catch(...) {} } { itk::NiftiImageIO::Pointer io = itk::NiftiImageIO::New(); if ( io->CanReadFile( path.c_str() ) ) { io->SetFileName( path.c_str() ); io->ReadImageInformation(); if ( io->GetPixelType() == itk::ImageIOBase::SCALAR && io->GetNumberOfDimensions()==4) { switch (io->GetDimensions(3)) { case 6: return true; break; case 15: return true; break; case 28: return true; break; case 45: return true; break; case 66: return true; break; case 91: return true; break; default : return false; } } } } return false; } DiffusionCoreIOMimeTypes::SHImageMimeType* DiffusionCoreIOMimeTypes::SHImageMimeType::Clone() const { return new SHImageMimeType(*this); } DiffusionCoreIOMimeTypes::SHImageMimeType DiffusionCoreIOMimeTypes::SH_MIMETYPE() { return SHImageMimeType(); } CustomMimeType DiffusionCoreIOMimeTypes::DTI_MIMETYPE() { CustomMimeType mimeType(DTI_MIMETYPE_NAME()); std::string category = "Tensor Image"; mimeType.SetComment("Diffusion Tensor Image"); mimeType.SetCategory(category); mimeType.AddExtension("dti"); return mimeType; } CustomMimeType DiffusionCoreIOMimeTypes::ODF_MIMETYPE() { CustomMimeType mimeType(ODF_MIMETYPE_NAME()); std::string category = "ODF Image"; mimeType.SetComment("Diffusion ODF Image"); mimeType.SetCategory(category); mimeType.AddExtension("odf"); mimeType.AddExtension("qbi"); // legacy support return mimeType; } // Names std::string DiffusionCoreIOMimeTypes::PEAK_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + "_PEAKS"; return name; } 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::DWI_DICOM_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".IMA"; return name; } std::string DiffusionCoreIOMimeTypes::DTI_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".dti"; return name; } std::string DiffusionCoreIOMimeTypes::ODF_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + ".odf"; return name; } std::string DiffusionCoreIOMimeTypes::SH_MIMETYPE_NAME() { static std::string name = IOMimeTypes::DEFAULT_BASE_NAME() + "_SH"; return name; } // Descriptions std::string DiffusionCoreIOMimeTypes::PEAK_MIMETYPE_DESCRIPTION() { static std::string description = "Peak Image"; return description; } 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::DWI_DICOM_MIMETYPE_DESCRIPTION() { static std::string description = "Diffusion Weighted Images"; return description; } std::string DiffusionCoreIOMimeTypes::DTI_MIMETYPE_DESCRIPTION() { static std::string description = "Diffusion Tensor Image"; return description; } std::string DiffusionCoreIOMimeTypes::ODF_MIMETYPE_DESCRIPTION() { static std::string description = "ODF Image"; return description; } std::string DiffusionCoreIOMimeTypes::SH_MIMETYPE_DESCRIPTION() { static std::string description = "SH Image"; return description; } } diff --git a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageReader.cpp b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageReader.cpp index 1fd94331b0..e0fe903cef 100644 --- a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageReader.cpp @@ -1,76 +1,83 @@ /*=================================================================== 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 "mitkPeakImageReader.h" #include #include #include #include #include #include #include #include #include #include #include namespace mitk { PeakImageReader::PeakImageReader(const PeakImageReader& other) : mitk::AbstractFileReader(other) { } PeakImageReader::PeakImageReader() : mitk::AbstractFileReader( CustomMimeType( mitk::DiffusionCoreIOMimeTypes::PEAK_MIMETYPE() ), mitk::DiffusionCoreIOMimeTypes::PEAK_MIMETYPE_DESCRIPTION() ) { m_ServiceReg = this->RegisterService(); } PeakImageReader::~PeakImageReader() { } std::vector > PeakImageReader::Read() { std::vector > result; + std::string location = GetInputLocation(); + std::string ext = itksys::SystemTools::GetFilenameExtension(location); typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(GetInputLocation()); + if (ext==".peak") + { + itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); + reader->SetImageIO(io); + } reader->Update(); Image::Pointer resultImage = dynamic_cast(PeakImage::New().GetPointer()); mitk::CastToMitkImage(reader->GetOutput(), resultImage); resultImage->SetVolume(reader->GetOutput()->GetBufferPointer()); StringProperty::Pointer nameProp; nameProp = StringProperty::New(itksys::SystemTools::GetFilenameWithoutExtension(GetInputLocation())); resultImage->SetProperty("name", nameProp); dynamic_cast(resultImage.GetPointer())->ConstructPolydata(); result.push_back( resultImage.GetPointer() ); return result; } } //namespace MITK mitk::PeakImageReader* mitk::PeakImageReader::Clone() const { return new PeakImageReader(*this); } diff --git a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageSerializer.cpp b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageSerializer.cpp index 87504a1ab5..b2f7168b15 100644 --- a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageSerializer.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkPeakImageSerializer.cpp @@ -1,80 +1,84 @@ /*=================================================================== 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 "mitkPeakImageSerializer.h" #include "mitkPeakImage.h" #include "itkImageFileWriter.h" #include #include - +#include MITK_REGISTER_SERIALIZER(PeakImageSerializer) mitk::PeakImageSerializer::PeakImageSerializer() { } mitk::PeakImageSerializer::~PeakImageSerializer() { } std::string mitk::PeakImageSerializer::Serialize() { const PeakImage* image = dynamic_cast( m_Data.GetPointer() ); if (image == nullptr) { MITK_ERROR << " Object at " << (const void*) this->m_Data << " is not an mitk::PeakImage. Cannot serialize as Nrrd."; return ""; } std::string filename( this->GetUniqueFilenameInWorkingDirectory() ); filename += "_"; filename += m_FilenameHint; - filename += ".nrrd"; + filename += ".peak"; std::string fullname(m_WorkingDirectory); fullname += "/"; fullname += itksys::SystemTools::ConvertToOutputPath(filename.c_str()); try { typedef mitk::ImageToItk< PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer itk_image = caster->GetOutput(); itk::ImageFileWriter< PeakImage::ItkPeakImageType >::Pointer writer = itk::ImageFileWriter< PeakImage::ItkPeakImageType >::New(); writer->SetInput(itk_image); writer->SetFileName(fullname); + itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); + io->SetFileType( itk::ImageIOBase::Binary ); + io->UseCompressionOn(); + writer->SetImageIO(io); writer->Update(); } catch (std::exception& e) { MITK_ERROR << " Error serializing object at " << (const void*) this->m_Data << " to " << fullname << ": " << e.what(); return ""; } return filename; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 8b5dbf5ee1..4d3adb37d6 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1023 +1,1024 @@ /*=================================================================== 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 "itkStreamlineTrackingFilter.h" #include #include #include #include #include "itkPointShell.h" #include #include #include #include #include #include #include #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_BuildFibersFinished(false) , m_BuildFibersReady(0) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_StoppingRegions(nullptr) , m_TargetRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_ExclusionRegions(nullptr) , m_OutputProbabilityMap(nullptr) , m_MinVoxelSize(-1) , m_AngularThresholdDeg(-1) , m_StepSizeVox(-1) , m_SamplingDistanceVox(-1) , m_AngularThreshold(-1) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) , m_AvoidStop(true) , m_RandomSampling(false) , m_SamplingDistance(-1) , m_DeflectionMod(1.0) , m_OnlyForwardSamples(true) , m_UseStopVotes(true) , m_NumberOfSamples(30) , m_NumPreviousDirections(1) , m_MaxNumTracts(-1) , m_Verbose(true) , m_LoopCheck(-1) , m_DemoMode(false) , m_Random(true) , m_UseOutputProbabilityMap(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) , m_InterpolateMasks(true) , m_TrialsPerSeed(10) , m_EndpointConstraint(EndpointConstraints::NONE) , m_IntroduceDirectionsFromPrior(true) , m_TrackingPriorAsMask(true) , m_TrackingPriorWeight(1.0) , m_TrackingPriorHandler(nullptr) { this->SetNumberOfRequiredInputs(0); } std::string StreamlineTrackingFilter::GetStatusText() { std::string status = "Seedpoints processed: " + boost::lexical_cast(m_Progress) + "/" + boost::lexical_cast(m_SeedPoints.size()); if (m_SeedPoints.size()>0) status += " (" + boost::lexical_cast(100*m_Progress/m_SeedPoints.size()) + "%)"; if (m_MaxNumTracts>0) status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts) + "/" + boost::lexical_cast(m_MaxNumTracts); else status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts); return status; } void StreamlineTrackingFilter::BeforeTracking() { m_StopTracking = false; m_TrackingHandler->SetRandom(m_Random); m_TrackingHandler->InitForTracking(); m_FiberPolyData = PolyDataType::New(); m_Points = vtkSmartPointer< vtkPoints >::New(); m_Cells = vtkSmartPointer< vtkCellArray >::New(); itk::Vector< double, 3 > imageSpacing = m_TrackingHandler->GetSpacing(); if(imageSpacing[0]SetAngularThreshold(m_AngularThreshold); if (m_TrackingPriorHandler!=nullptr) { m_TrackingPriorHandler->SetRandom(m_Random); m_TrackingPriorHandler->InitForTracking(); m_TrackingPriorHandler->SetAngularThreshold(m_AngularThreshold); } if (m_SamplingDistanceVoxGetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } if (m_UseOutputProbabilityMap) { m_OutputProbabilityMap = ItkDoubleImgType::New(); m_OutputProbabilityMap->SetSpacing(imageSpacing); m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin()); m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection()); m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion()); m_OutputProbabilityMap->Allocate(); m_OutputProbabilityMap->FillBuffer(0); } m_MaskInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_StopInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_SeedInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_TargetInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_ExclusionInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); if (m_StoppingRegions.IsNull()) { m_StoppingRegions = ItkFloatImgType::New(); m_StoppingRegions->SetSpacing( imageSpacing ); m_StoppingRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_StoppingRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_StoppingRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_StoppingRegions->Allocate(); m_StoppingRegions->FillBuffer(0); } else std::cout << "StreamlineTracking - Using stopping region image" << std::endl; m_StopInterpolator->SetInputImage(m_StoppingRegions); if (m_ExclusionRegions.IsNotNull()) { std::cout << "StreamlineTracking - Using exclusion region image" << std::endl; m_ExclusionInterpolator->SetInputImage(m_ExclusionRegions); } if (m_TargetRegions.IsNull()) { m_TargetImageSet = false; m_TargetRegions = ItkFloatImgType::New(); m_TargetRegions->SetSpacing( imageSpacing ); m_TargetRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_TargetRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_TargetRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_TargetRegions->Allocate(); m_TargetRegions->FillBuffer(1); } else { m_TargetImageSet = true; m_TargetInterpolator->SetInputImage(m_TargetRegions); std::cout << "StreamlineTracking - Using target region image" << std::endl; } if (m_SeedImage.IsNull()) { m_SeedImageSet = false; m_SeedImage = ItkFloatImgType::New(); m_SeedImage->SetSpacing( imageSpacing ); m_SeedImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_SeedImage->SetDirection( m_TrackingHandler->GetDirection() ); m_SeedImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_SeedImage->Allocate(); m_SeedImage->FillBuffer(1); } else { m_SeedImageSet = true; std::cout << "StreamlineTracking - Using seed image" << std::endl; } m_SeedInterpolator->SetInputImage(m_SeedImage); if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkFloatImgType::New(); m_MaskImage->SetSpacing( imageSpacing ); m_MaskImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_MaskImage->SetDirection( m_TrackingHandler->GetDirection() ); m_MaskImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } else std::cout << "StreamlineTracking - Using mask image" << std::endl; m_MaskInterpolator->SetInputImage(m_MaskImage); // Autosettings for endpoint constraints if (m_EndpointConstraint==EndpointConstraints::NONE && m_TargetImageSet && m_SeedImageSet) { MITK_INFO << "No endpoint constraint chosen but seed and target image set --> setting constraint to EPS_IN_SEED_AND_TARGET"; m_EndpointConstraint = EndpointConstraints::EPS_IN_SEED_AND_TARGET; } else if (m_EndpointConstraint==EndpointConstraints::NONE && m_TargetImageSet) { MITK_INFO << "No endpoint constraint chosen but target image set --> setting constraint to EPS_IN_TARGET"; m_EndpointConstraint = EndpointConstraints::EPS_IN_TARGET; } // Check if endpoint constraints are valid FiberType test_fib; itk::Point p; p.Fill(0); test_fib.push_back(p); test_fib.push_back(p); IsValidFiber(&test_fib); if (m_SeedPoints.empty()) GetSeedPointsFromSeedImage(); m_BuildFibersReady = 0; m_BuildFibersFinished = false; m_Tractogram.clear(); m_SamplingPointset = mitk::PointSet::New(); m_AlternativePointset = mitk::PointSet::New(); m_StopVotePointset = mitk::PointSet::New(); m_StartTime = std::chrono::system_clock::now(); if (m_DemoMode) omp_set_num_threads(1); if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::DETERMINISTIC) std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; else if(m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::PROBABILISTIC) { std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; std::cout << "StreamlineTracking - Trials per seed: " << m_TrialsPerSeed << std::endl; } else std::cout << "StreamlineTracking - Mode: ???" << std::endl; if (m_EndpointConstraint==EndpointConstraints::NONE) std::cout << "StreamlineTracking - Endpoint constraint: NONE" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET_LABELDIFF" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_SEED_AND_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_SEED_AND_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::MIN_ONE_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: MIN_ONE_EP_IN_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::ONE_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: ONE_EP_IN_TARGET" << std::endl; else if (m_EndpointConstraint==EndpointConstraints::NO_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: NO_EP_IN_TARGET" << std::endl; std::cout << "StreamlineTracking - Angular threshold: " << m_AngularThreshold << " (" << 180*std::acos( m_AngularThreshold )/itk::Math::pi << "°)" << std::endl; std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/m_MinVoxelSize << "*vox)" << std::endl; std::cout << "StreamlineTracking - Seeds per voxel: " << m_SeedsPerVoxel << std::endl; std::cout << "StreamlineTracking - Max. tract length: " << m_MaxTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Min. tract length: " << m_MinTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Max. num. tracts: " << m_MaxNumTracts << std::endl; std::cout << "StreamlineTracking - Loop check: " << m_LoopCheck << "°" << std::endl; std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_NumberOfSamples << std::endl; std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/m_MinVoxelSize << "*vox)" << std::endl; std::cout << "StreamlineTracking - Deflection modifier: " << m_DeflectionMod << std::endl; std::cout << "StreamlineTracking - Use stop votes: " << m_UseStopVotes << std::endl; std::cout << "StreamlineTracking - Only frontal samples: " << m_OnlyForwardSamples << std::endl; if (m_TrackingPriorHandler!=nullptr) - std::cout << "StreamlineTracking - Using directional prior for tractography" << std::endl; + std::cout << "StreamlineTracking - Using directional prior for tractography (w=" << m_TrackingPriorWeight << ")" << std::endl; if (m_DemoMode) { std::cout << "StreamlineTracking - Running in demo mode"; std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl; } else std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl; } void StreamlineTrackingFilter::CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir) { pos[0] += dir[0]*m_StepSize; pos[1] += dir[1]*m_StepSize; pos[2] += dir[2]*m_StepSize; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< float > theta; theta.resize(NPoints); std::vector< float > phi; phi.resize(NPoints); float C = sqrt(4*itk::Math::pi); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(int i=0; i0 && i d; d[0] = cos(theta[i]) * cos(phi[i]); d[1] = cos(theta[i]) * sin(phi[i]); d[2] = sin(theta[i]); pointshell.push_back(d); } return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(const itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { if (m_DemoMode) { m_SamplingPointset->Clear(); m_AlternativePointset->Clear(); m_StopVotePointset->Clear(); } vnl_vector_fixed direction; direction.fill(0); if (mitk::imv::IsInsideMask(pos, m_InterpolateMasks, m_MaskInterpolator) && !mitk::imv::IsInsideMask(pos, m_InterpolateMasks, m_StopInterpolator)) direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position else return direction; int stop_votes = 0; int possible_stop_votes = 0; if (!olddirs.empty()) { vnl_vector_fixed olddir = olddirs.back(); std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_NumberOfSamples); itk::Point sample_pos; int alternatives = 1; for (unsigned int i=0; i d; bool is_stop_voter = false; if (m_Random && m_RandomSampling) { d[0] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[1] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[2] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d.normalize(); d *= m_TrackingHandler->GetRandDouble(0,m_SamplingDistance); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); if (m_UseStopVotes && dot>0.7) { is_stop_voter = true; possible_stop_votes++; } else if (m_OnlyForwardSamples && dot<0) continue; d *= m_SamplingDistance; } sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (mitk::imv::IsInsideMask(sample_pos, m_InterpolateMasks, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) { direction += tempDir; if(m_DemoMode) m_SamplingPointset->InsertPoint(i, sample_pos); } else if (m_AvoidStop && olddir.magnitude()>0.5) // out of white matter { if (is_stop_voter) stop_votes++; if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); float dot = dot_product(d, olddir); if (dot >= 0.0) // in front of plane defined by pos and olddir d = -d + 2*dot*olddir; // reflect else d = -d; // invert // look a bit further into the other direction sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; alternatives++; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (mitk::imv::IsInsideMask(sample_pos, m_InterpolateMasks, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) // are we back in the white matter? { direction += d * m_DeflectionMod; // go into the direction of the white matter direction += tempDir; // go into the direction of the white matter direction at this location if(m_DemoMode) m_AlternativePointset->InsertPoint(alternatives, sample_pos); } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); } } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); if (is_stop_voter) stop_votes++; } } } bool valid = false; if (direction.magnitude()>0.001 && (possible_stop_votes==0 || (float)stop_votes/possible_stop_votes<0.5) ) { direction.normalize(); valid = true; } else direction.fill(0); if (m_TrackingPriorHandler!=nullptr && (m_IntroduceDirectionsFromPrior || valid)) { vnl_vector_fixed prior = m_TrackingPriorHandler->ProposeDirection(pos, olddirs, oldIndex); if (prior.magnitude()>0.001) { prior.normalize(); if (dot_product(prior,direction)<0) prior *= -1; direction = (1.0f-m_TrackingPriorWeight) * direction + m_TrackingPriorWeight * prior; direction.normalize(); } else if (m_TrackingPriorAsMask) direction.fill(0.0); } return direction; } float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, DirectionContainer* container, float tractLength, bool front, bool &exclude) { vnl_vector_fixed zero_dir; zero_dir.fill(0.0); std::deque< vnl_vector_fixed > last_dirs; for (unsigned int i=0; i oldIndex; m_TrackingHandler->WorldToIndex(pos, oldIndex); // get new position CalculateNewPosition(pos, dir); if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask(pos, m_InterpolateMasks, m_ExclusionInterpolator)) { exclude = true; return tractLength; } if (m_AbortTracking) return tractLength; // if yes, add new point to streamline dir.normalize(); if (front) { fib->push_front(pos); container->push_front(dir); } else { fib->push_back(pos); container->push_back(dir); } tractLength += m_StepSize; if (m_LoopCheck>=0 && CheckCurvature(container, front)>m_LoopCheck) return tractLength; if (tractLength>m_MaxTractLength) return tractLength; if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? { #pragma omp critical { m_BuildFibersReady++; m_Tractogram.push_back(*fib); BuildFibers(true); m_Stop = true; while (m_Stop){ } } } last_dirs.push_back(dir); if (last_dirs.size()>m_NumPreviousDirections) last_dirs.pop_front(); dir = GetNewDirection(pos, last_dirs, oldIndex); while (m_PauseTracking){} if (dir.magnitude()<0.0001) return tractLength; } return tractLength; } float StreamlineTrackingFilter::CheckCurvature(DirectionContainer* fib, bool front) { if (fib->size()<8) return 0; float m_Distance = std::max(m_MinVoxelSize*4, m_StepSize*8); float dist = 0; std::vector< vnl_vector_fixed< float, 3 > > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); float dev = 0; if (front) { int c = 0; while(distsize()-1) { dist += m_StepSize; vnl_vector_fixed< float, 3 > v = fib->at(c); + if (dot_product(v,meanV)<0) + v = -v; vectors.push_back(v); meanV += v; c++; } } else { int c = fib->size()-1; while(dist=0) { dist += m_StepSize; vnl_vector_fixed< float, 3 > v = fib->at(c); + if (dot_product(v,meanV)<0) + v = -v; vectors.push_back(v); meanV += v; c--; } } meanV.normalize(); for (unsigned int c=0; c1.0) angle = 1.0; dev += acos(angle)*180/itk::Math::pi; } if (vectors.size()>0) dev /= vectors.size(); return dev; } void StreamlineTrackingFilter::SetTrackingPriorHandler(mitk::TrackingDataHandler *TrackingPriorHandler) { m_TrackingPriorHandler = TrackingPriorHandler; } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - Calculating seed points."; m_SeedPoints.clear(); typedef ImageRegionConstIterator< ItkFloatImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); sit.GoToBegin(); while (!sit.IsAtEnd()) { if (sit.Value()>0) { ItkFloatImgType::IndexType index = sit.GetIndex(); itk::ContinuousIndex start; start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); if ( mitk::imv::IsInsideMask(worldPos, m_InterpolateMasks, m_MaskInterpolator) ) { m_SeedPoints.push_back(worldPos); for (int s = 1; s < m_SeedsPerVoxel; s++) { start[0] = index[0] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[1] = index[1] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[2] = index[2] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); m_SeedPoints.push_back(worldPos); } } } ++sit; } } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (m_Random) std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); m_CurrentTracts = 0; int num_seeds = m_SeedPoints.size(); itk::Index<3> zeroIndex; zeroIndex.Fill(0); m_Progress = 0; int i = 0; int print_interval = num_seeds/100; if (print_interval<100) m_Verbose=false; #pragma omp parallel while (i=num_seeds || m_StopTracking) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { m_Progress += print_interval; std::cout << " \r"; if (m_MaxNumTracts>0) std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_MaxNumTracts << '\r'; else std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << '\r'; cout.flush(); } const itk::Point worldPos = m_SeedPoints.at(temp_i); for (unsigned int trials=0; trials dir; dir.fill(0.0); std::deque< vnl_vector_fixed > olddirs; dir = GetNewDirection(worldPos, olddirs, zeroIndex) * 0.5f; bool exclude = false; if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask(worldPos, m_InterpolateMasks, m_ExclusionInterpolator)) exclude = true; bool success = false; if (dir.magnitude()>0.0001 && !exclude) { // forward tracking tractLength = FollowStreamline(worldPos, dir, &fib, &direction_container, 0, false, exclude); fib.push_front(worldPos); // backward tracking if (!exclude) tractLength = FollowStreamline(worldPos, -dir, &fib, &direction_container, tractLength, true, exclude); counter = fib.size(); if (tractLength>=m_MinTractLength && counter>=2 && !exclude) { #pragma omp critical if ( IsValidFiber(&fib) ) { if (!m_StopTracking) { if (!m_UseOutputProbabilityMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); m_CurrentTracts++; success = true; } if (m_MaxNumTracts > 0 && m_CurrentTracts>=static_cast(m_MaxNumTracts)) { if (!m_StopTracking) { std::cout << " \r"; MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; } m_StopTracking = true; } } } } if (success || m_TrackingHandler->GetMode()!=mitk::TrackingDataHandler::PROBABILISTIC) break; // we only try one seed point multiple times if we use a probabilistic tracker and have not found a valid streamline yet }// trials per seed }// seed points this->AfterTracking(); } bool StreamlineTrackingFilter::IsValidFiber(FiberType* fib) { if (m_EndpointConstraint==EndpointConstraints::NONE) { return true; } else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) { if (m_TargetImageSet) { float v1 = mitk::imv::GetImageValue(fib->front(), false, m_TargetInterpolator); float v2 = mitk::imv::GetImageValue(fib->back(), false, m_TargetInterpolator); if ( v1>0.0 && v2>0.0 && v1!=v2 ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET_LABELDIFF chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_SEED_AND_TARGET) { if (m_TargetImageSet && m_SeedImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; if ( mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target or seed image set but endpoint constraint EPS_IN_SEED_AND_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::MIN_ONE_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) || mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint MIN_ONE_EP_IN_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::ONE_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) && !mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return true; if ( !mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) && mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint ONE_EP_IN_TARGET chosen!"; } else if (m_EndpointConstraint==EndpointConstraints::NO_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) || mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) return false; return true; } else mitkThrow() << "No target image set but endpoint constraint NO_EP_IN_TARGET chosen!"; } return true; } void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib) { ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0); for (auto p : *fib) { ItkDoubleImgType::IndexType idx; m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); if (idx != last_idx) { if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(idx)) m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1); last_idx = idx; } } } void StreamlineTrackingFilter::BuildFibers(bool check) { if (m_BuildFibersReady::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); for (unsigned int i=0; i container = vtkSmartPointer::New(); FiberType fib = m_Tractogram.at(i); for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) { vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if (check) for (int i=0; iSetPoints(vNewPoints); m_FiberPolyData->SetLines(vNewLines); m_BuildFibersFinished = true; } void StreamlineTrackingFilter::AfterTracking() { if (m_Verbose) std::cout << " \r"; if (!m_UseOutputProbabilityMap) { MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; MITK_INFO << "Generating polydata "; BuildFibers(false); } else { itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer filter = itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); filter->SetInput(m_OutputProbabilityMap); filter->SetOutputMaximum(1.0); filter->SetOutputMinimum(0.0); filter->Update(); m_OutputProbabilityMap = filter->GetOutput(); } MITK_INFO << "done"; m_EndTime = std::chrono::system_clock::now(); std::chrono::hours hh = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::minutes mm = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::seconds ss = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; ss %= 60; MITK_INFO << "Tracking took " << hh.count() << "h, " << mm.count() << "m and " << ss.count() << "s"; m_SeedPoints.clear(); - - if (m_TrackingPriorHandler!=nullptr) - delete m_TrackingPriorHandler; } void StreamlineTrackingFilter::SetDicomProperties(mitk::FiberBundle::Pointer fib) { std::string model_code_value = "-"; std::string model_code_meaning = "-"; std::string algo_code_value = "-"; std::string algo_code_meaning = "-"; if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC && dynamic_cast(m_TrackingHandler) && !m_TrackingHandler->GetInterpolate()) { algo_code_value = "sup181_ee04"; algo_code_meaning = "FACT"; } else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC) { algo_code_value = "sup181_ee01"; algo_code_meaning = "Deterministic"; } else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::PROBABILISTIC) { algo_code_value = "sup181_ee02"; algo_code_meaning = "Probabilistic"; } if (dynamic_cast(m_TrackingHandler) || (dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetIsOdfFromTensor() ) ) { if ( dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetNumTensorImages()>1 ) { model_code_value = "sup181_bb02"; model_code_meaning = "Multi Tensor"; } else { model_code_value = "sup181_bb01"; model_code_meaning = "Single Tensor"; } } else if (dynamic_cast*>(m_TrackingHandler) || dynamic_cast*>(m_TrackingHandler)) { model_code_value = "sup181_bb03"; model_code_meaning = "Model Free"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "ODF"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "Peaks"; } fib->SetProperty("DICOM.anatomy.value", mitk::StringProperty::New("T-A0095")); fib->SetProperty("DICOM.anatomy.meaning", mitk::StringProperty::New("White matter of brain and spinal cord")); fib->SetProperty("DICOM.algo_code.value", mitk::StringProperty::New(algo_code_value)); fib->SetProperty("DICOM.algo_code.meaning", mitk::StringProperty::New(algo_code_meaning)); fib->SetProperty("DICOM.model_code.value", mitk::StringProperty::New(model_code_value)); fib->SetProperty("DICOM.model_code.meaning", mitk::StringProperty::New(model_code_meaning)); } } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp index 9968c4979c..425e5e764b 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp @@ -1,520 +1,574 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include +#include #define _USE_MATH_DEFINES #include const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("input", "i", mitkCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false); parser.addArgument("algorithm", "a", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle/probability map", us::Any(), false); parser.endGroup(); parser.beginGroup("2. Seeding:"); parser.addArgument("seeds", "", mitkCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_image", "", mitkCommandLineParser::String, "Seed image:", "mask image defining seed voxels", us::Any()); parser.addArgument("trials_per_seed", "", mitkCommandLineParser::Int, "Max. trials per seed:", "try each seed N times until a valid streamline is obtained (only for probabilistic tractography)", 10); parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded", -1); parser.endGroup(); parser.beginGroup("3. Tractography constraints:"); parser.addArgument("tracking_mask", "", mitkCommandLineParser::String, "Mask image:", "streamlines leaving the mask will stop immediately", us::Any()); parser.addArgument("stop_image", "", mitkCommandLineParser::String, "Stop ROI image:", "streamlines entering the mask will stop immediately", us::Any()); parser.addArgument("exclusion_image", "", mitkCommandLineParser::String, "Exclusion ROI image:", "streamlines entering the mask will be discarded", us::Any()); parser.addArgument("ep_constraint", "", mitkCommandLineParser::String, "Endpoint constraint:", "determines which fibers are accepted based on their endpoint location - options are NONE, EPS_IN_TARGET, EPS_IN_TARGET_LABELDIFF, EPS_IN_SEED_AND_TARGET, MIN_ONE_EP_IN_TARGET, ONE_EP_IN_TARGET and NO_EP_IN_TARGET", us::Any()); parser.addArgument("target_image", "", mitkCommandLineParser::String, "Target ROI image:", "effact depends on the chosen endpoint constraint (option ep_constraint)", us::Any()); parser.endGroup(); parser.beginGroup("4. Streamline integration parameters:"); parser.addArgument("sharpen_odfs", "", mitkCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkCommandLineParser::Float, "ODF Cutoff:", "threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.0); parser.addArgument("step_size", "", mitkCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("min_tract_length", "", mitkCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("angular_threshold", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size, minimum 15°)"); parser.addArgument("loop_check", "", mitkCommandLineParser::Float, "Check for loops:", "threshold on angular stdev over the last 4 voxel lengths"); parser.endGroup(); - parser.beginGroup("5. Neighborhood sampling:"); + parser.beginGroup("5. Tractography prior:"); + parser.addArgument("prior_image", "", mitkCommandLineParser::String, "Peak prior:", "tractography prior in thr for of a peak image", us::Any()); + parser.addArgument("prior_weight", "", mitkCommandLineParser::Float, "Prior weight", "weighting factor between prior and data.", 0.5); + parser.addArgument("restrict_to_prior", "", mitkCommandLineParser::Bool, "Restrict to prior:", "restrict tractography to regions where the prior is valid."); + parser.addArgument("new_directions_from_prior", "", mitkCommandLineParser::Bool, "New directios from prior:", "the prior can create directions where there are none in the data."); + parser.endGroup(); + + parser.beginGroup("6. Neighborhood sampling:"); parser.addArgument("num_samples", "", mitkCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.endGroup(); - parser.beginGroup("6. Tensor tractography specific:"); + parser.beginGroup("7. Tensor tractography specific:"); parser.addArgument("tend_f", "", mitkCommandLineParser::Float, "Weight f", "weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkCommandLineParser::Float, "Weight g", "weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); parser.endGroup(); - parser.beginGroup("7. Random forest tractography specific:"); + parser.beginGroup("8. Random forest tractography specific:"); parser.addArgument("forest", "", mitkCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any()); parser.addArgument("use_sh_features", "", mitkCommandLineParser::Bool, "Use SH features:", "use SH features"); parser.endGroup(); - parser.beginGroup("8. Additional input:"); + parser.beginGroup("9. Additional input:"); parser.addArgument("additional_images", "", mitkCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional features for RF tractography)", us::Any()); parser.endGroup(); - parser.beginGroup("9. Misc:"); + parser.beginGroup("10. Misc:"); parser.addArgument("flip_x", "", mitkCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); parser.addArgument("no_data_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate input data:", "don't interpolate input image values"); parser.addArgument("no_mask_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate masks:", "don't interpolate mask image values"); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "compress output fibers using the given error threshold (in mm)"); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["input"]); std::string outFile = us::any_cast(parsedArgs["out"]); std::string algorithm = us::any_cast(parsedArgs["algorithm"]); + std::string prior_image = ""; + if (parsedArgs.count("prior_image")) + prior_image = us::any_cast(parsedArgs["prior_image"]); + + float prior_weight = 0.5; + if (parsedArgs.count("prior_weight")) + prior_weight = us::any_cast(parsedArgs["prior_weight"]); + + bool restrict_to_prior = false; + if (parsedArgs.count("restrict_to_prior")) + restrict_to_prior = us::any_cast(parsedArgs["restrict_to_prior"]); + + bool new_directions_from_prior = false; + if (parsedArgs.count("new_directions_from_prior")) + new_directions_from_prior = us::any_cast(parsedArgs["new_directions_from_prior"]); + bool sharpen_odfs = false; if (parsedArgs.count("sharpen_odfs")) sharpen_odfs = us::any_cast(parsedArgs["sharpen_odfs"]); bool interpolate = true; if (parsedArgs.count("no_data_interpolation")) interpolate = !us::any_cast(parsedArgs["no_data_interpolation"]); bool mask_interpolation = true; if (parsedArgs.count("no_mask_interpolation")) interpolate = !us::any_cast(parsedArgs["no_mask_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); bool use_stop_votes = false; if (parsedArgs.count("use_stop_votes")) use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); bool use_only_forward_samples = false; if (parsedArgs.count("use_only_forward_samples")) use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); bool flip_x = false; if (parsedArgs.count("flip_x")) flip_x = us::any_cast(parsedArgs["flip_x"]); bool flip_y = false; if (parsedArgs.count("flip_y")) flip_y = us::any_cast(parsedArgs["flip_y"]); bool flip_z = false; if (parsedArgs.count("flip_z")) flip_z = us::any_cast(parsedArgs["flip_z"]); bool apply_image_rotation = false; if (parsedArgs.count("apply_image_rotation")) apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); float min_tract_length = 20; if (parsedArgs.count("min_tract_length")) min_tract_length = us::any_cast(parsedArgs["min_tract_length"]); float loop_check = -1; if (parsedArgs.count("loop_check")) loop_check = us::any_cast(parsedArgs["loop_check"]); std::string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); std::string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); std::string seedFile = ""; if (parsedArgs.count("seed_image")) seedFile = us::any_cast(parsedArgs["seed_image"]); std::string targetFile = ""; if (parsedArgs.count("target_image")) targetFile = us::any_cast(parsedArgs["target_image"]); std::string exclusionFile = ""; if (parsedArgs.count("exclusion_image")) exclusionFile = us::any_cast(parsedArgs["exclusion_image"]); std::string stopFile = ""; if (parsedArgs.count("stop_image")) stopFile = us::any_cast(parsedArgs["stop_image"]); std::string ep_constraint = "NONE"; if (parsedArgs.count("ep_constraint")) ep_constraint = us::any_cast(parsedArgs["ep_constraint"]); float cutoff = 0.1; if (parsedArgs.count("cutoff")) cutoff = us::any_cast(parsedArgs["cutoff"]); float odf_cutoff = 0.0; if (parsedArgs.count("odf_cutoff")) odf_cutoff = us::any_cast(parsedArgs["odf_cutoff"]); float stepsize = -1; if (parsedArgs.count("step_size")) stepsize = us::any_cast(parsedArgs["step_size"]); float sampling_distance = -1; if (parsedArgs.count("sampling_distance")) sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); int num_samples = 0; if (parsedArgs.count("num_samples")) num_samples = us::any_cast(parsedArgs["num_samples"]); int num_seeds = 1; if (parsedArgs.count("seeds")) num_seeds = us::any_cast(parsedArgs["seeds"]); unsigned int trials_per_seed = 10; if (parsedArgs.count("trials_per_seed")) trials_per_seed = us::any_cast(parsedArgs["trials_per_seed"]); float tend_f = 1; if (parsedArgs.count("tend_f")) tend_f = us::any_cast(parsedArgs["tend_f"]); float tend_g = 0; if (parsedArgs.count("tend_g")) tend_g = us::any_cast(parsedArgs["tend_g"]); float angular_threshold = -1; if (parsedArgs.count("angular_threshold")) angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); unsigned int max_tracts = -1; if (parsedArgs.count("max_tracts")) max_tracts = us::any_cast(parsedArgs["max_tracts"]); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (ext != ".fib" && ext != ".trk") { MITK_INFO << "Output file format not supported. Use one of .fib, .trk, .nii, .nii.gz, .nrrd"; return EXIT_FAILURE; } // LOAD DATASETS mitkCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkFloatImgType; MITK_INFO << "loading input"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(mitk::IOUtil::Load(input_files.at(i))[0].GetPointer()); input_images.push_back(mitkImage); } ItkFloatImgType::Pointer mask = nullptr; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(maskFile)[0].GetPointer()); mask = ItkFloatImgType::New(); mitk::CastToItkImage(img, mask); } ItkFloatImgType::Pointer seed = nullptr; if (!seedFile.empty()) { MITK_INFO << "loading seed ROI image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(seedFile)[0].GetPointer()); seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); } ItkFloatImgType::Pointer stop = nullptr; if (!stopFile.empty()) { MITK_INFO << "loading stop ROI image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(stopFile)[0].GetPointer()); stop = ItkFloatImgType::New(); mitk::CastToItkImage(img, stop); } ItkFloatImgType::Pointer target = nullptr; if (!targetFile.empty()) { MITK_INFO << "loading target ROI image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(targetFile)[0].GetPointer()); target = ItkFloatImgType::New(); mitk::CastToItkImage(img, target); } ItkFloatImgType::Pointer exclusion = nullptr; if (!exclusionFile.empty()) { MITK_INFO << "loading exclusion ROI image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(exclusionFile)[0].GetPointer()); exclusion = ItkFloatImgType::New(); mitk::CastToItkImage(img, exclusion); } MITK_INFO << "loading additional images"; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(file)[0].GetPointer()); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); if (algorithm == "ProbTensor") { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkTensorImg = caster->GetOutput(); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkTensorImg ); filter->Update(); mitk::Image::Pointer image = mitk::Image::New(); FilterType::OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); input_images.clear(); input_images.push_back(image); sharpen_odfs = true; odf_cutoff = 0; } typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); + if (!prior_image.empty()) + { + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); + mitk::PeakImage::Pointer priorImage = dynamic_cast(mitk::IOUtil::Load(prior_image, &functor)[0].GetPointer()); + + if (priorImage.IsNull()) + { + MITK_INFO << "Only peak images are supported as prior at the moment!"; + return EXIT_FAILURE; + } + + mitk::TrackingDataHandler* priorhandler = new mitk::TrackingHandlerPeaks(); + + typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(priorImage); + caster->Update(); + mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); + + dynamic_cast(priorhandler)->SetPeakImage(itkImg); + dynamic_cast(priorhandler)->SetPeakThreshold(0.0); + dynamic_cast(priorhandler)->SetInterpolate(interpolate); + dynamic_cast(priorhandler)->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); + + tracker->SetTrackingPriorHandler(priorhandler); + tracker->SetTrackingPriorWeight(prior_weight); + tracker->SetTrackingPriorAsMask(restrict_to_prior); + tracker->SetIntroduceDirectionsFromPrior(new_directions_from_prior); + } + mitk::TrackingDataHandler* handler; if (algorithm == "DetRF" || algorithm == "ProbRF") { mitk::TractographyForest::Pointer forest = dynamic_cast(mitk::IOUtil::Load(forestFile)[0].GetPointer()); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } if (algorithm == "ProbRF") handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); } else if (algorithm == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetPeakImage(itkImg); dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); dynamic_cast(handler)->SetPeakThreshold(cutoff); } else if (algorithm == "DetTensor") { handler = new mitk::TrackingHandlerTensor(); for (auto input_image : input_images) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_image); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(handler)->AddTensorImage(itkImg); } dynamic_cast(handler)->SetFaThreshold(cutoff); dynamic_cast(handler)->SetF(tend_f); dynamic_cast(handler)->SetG(tend_g); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (algorithm == "DetODF" || algorithm == "ProbODF" || algorithm == "ProbTensor") { handler = new mitk::TrackingHandlerOdf(); typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetOdfImage(itkImg); dynamic_cast(handler)->SetGfaThreshold(cutoff); dynamic_cast(handler)->SetOdfThreshold(odf_cutoff); dynamic_cast(handler)->SetSharpenOdfs(sharpen_odfs); if (algorithm == "ProbODF" || algorithm == "ProbTensor") dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); if (algorithm == "ProbTensor") dynamic_cast(handler)->SetIsOdfFromTensor(true); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { MITK_INFO << "Unknown tractography algorithm (" + algorithm+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } handler->SetInterpolate(interpolate); handler->SetFlipX(flip_x); handler->SetFlipY(flip_y); handler->SetFlipZ(flip_z); if (ep_constraint=="NONE") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); else if (ep_constraint=="EPS_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); else if (ep_constraint=="EPS_IN_TARGET_LABELDIFF") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); else if (ep_constraint=="EPS_IN_SEED_AND_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); else if (ep_constraint=="MIN_ONE_EP_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); else if (ep_constraint=="ONE_EP_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); else if (ep_constraint=="NO_EP_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); MITK_INFO << "Tractography algorithm: " << algorithm; tracker->SetInterpolateMasks(mask_interpolation); tracker->SetNumberOfSamples(num_samples); tracker->SetAngularThreshold(angular_threshold); tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetTargetRegions(target); tracker->SetExclusionRegions(exclusion); tracker->SetSeedsPerVoxel(num_seeds); tracker->SetStepSize(stepsize); tracker->SetSamplingDistance(sampling_distance); tracker->SetUseStopVotes(use_stop_votes); tracker->SetOnlyForwardSamples(use_only_forward_samples); tracker->SetLoopCheck(loop_check); tracker->SetMaxNumTracts(max_tracts); tracker->SetTrialsPerSeed(trials_per_seed); tracker->SetTrackingHandler(handler); if (ext != ".fib" && ext != ".trk") tracker->SetUseOutputProbabilityMap(true); tracker->SetMinTractLength(min_tract_length); tracker->Update(); if (ext == ".fib" || ext == ".trk") { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui index 8435a5d6bd..539ed2cc15 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui @@ -1,395 +1,411 @@ QmitkFiberQuantificationViewControls 0 0 365 581 Form 25 - - + + - Input Data + Fiber-derived images - + 0 0 0 0 - - - - - - - - + + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Perform selected operation on all selected fiber bundles. + - Tractogram: + Generate Image - - - - Reference Image: + + + + + 0 + 0 + + + + Upsampling factor + + + 1 + + + 0.100000000000000 + + + 10.000000000000000 + + + 0.100000000000000 + + + 1.000000000000000 + + + + + 0 + 0 + + + + + Tract Density Image (TDI) + + + + + Normalized TDI + + + + + Binary Envelope + + + + + Fiber Bundle Image + + + + + Fiber Endings Image + + + + + Fiber Endings Pointset + + + + Principal Fiber Directions 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 Fiber directions with an angle smaller than the defined threshold are clustered. 2 0.000000000000000 90.000000000000000 1.000000000000000 30.000000000000000 0 0 <html><head/><body><p>Directions shorter than the defined threshold are discarded.</p></body></html> 3 1.000000000000000 0.100000000000000 0.300000000000000 Angular Threshold: Max. clusters: Size Threshold: 0 0 Maximum number of fiber directions per voxel. 100 3 Normalization: 0 0 + + 0 + Global maximum Single vector Voxel-wise maximum 0 0 Image containing the number of distinct fiber clusters per voxel. Output #Directions per Voxel - true + false false Generate Directions - - + + - Fiber-derived images + Input Data - + 0 0 0 0 - - - - false - - - - 0 - 0 - - - - - 200 - 16777215 - - - - - 11 - - - - Perform selected operation on all selected fiber bundles. - - - Generate Image - - - - - - - 0 - 0 - - - - Upsampling factor - - - 1 - - - 0.100000000000000 - - - 10.000000000000000 - - - 0.100000000000000 - - - 1.000000000000000 + + + + + + + + + Tractogram: - - - - - 0 - 0 - + + + + Reference Image: - - - Tract Density Image (TDI) - - - - - Normalized TDI - - - - - Binary Envelope - - - - - Fiber Bundle Image - - - - - Fiber Endings Image - - - - - Fiber Endings Pointset - - + + + + Qt::Vertical + + + + 20 + 40 + + + + QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index ba2f299036..4877541e23 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,969 +1,971 @@ /*=================================================================== 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 #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include // VTK #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) : m_View(view) { } void QmitkStreamlineTrackingWorker::run() { m_View->m_Tracker->Update(); m_View->m_TrackingThread.quit(); } QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_TrackingWorker(this) , m_Controls(nullptr) , m_FirstTensorProbRun(true) , m_FirstInteractiveRun(true) , m_TrackingHandler(nullptr) , m_ThreadIsRunning(false) , m_DeleteTrackingHandler(false) , m_Visible(false) + , m_LastPrior("") + , m_TrackingPriorHandler(nullptr) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TargetImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_PriorImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ForestBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ExclusionImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isPeakImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestBox->SetPredicate(isTractographyForest); m_Controls->m_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageBox->SetZeroEntryText("--"); m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_SeedImageBox->SetZeroEntryText("--"); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_StopImageBox->SetZeroEntryText("--"); m_Controls->m_TargetImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_TargetImageBox->SetZeroEntryText("--"); m_Controls->m_PriorImageBox->SetPredicate( isPeakImagePredicate ); m_Controls->m_PriorImageBox->SetZeroEntryText("--"); m_Controls->m_ExclusionImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_ExclusionImageBox->SetZeroEntryText("--"); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_SeedImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TargetImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_PriorImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ExclusionImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(valueChanged(double)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(valueChanged(int)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskInterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_LoopCheckBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TrialsPerSeedBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_EpConstraintsBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); StartStopTrackingGui(false); } UpdateGui(); } void QmitkStreamlineTrackingView::StopTractography() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); } void QmitkStreamlineTrackingView::TimerUpdate() { if (m_Tracker.IsNull()) return; QString status_text(m_Tracker->GetStatusText().c_str()); m_Controls->m_StatusTextBox->setText(status_text); } void QmitkStreamlineTrackingView::BeforeThread() { m_TrackingTimer->start(1000); } void QmitkStreamlineTrackingView::AfterThread() { m_TrackingTimer->stop(); if (!m_Tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_ParentNode->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked() && fiberBundle->GetNumberOfLines()>0) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); m_Tracker->SetDicomProperties(fib); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_ParentNode->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); node->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); GetDataStorage()->Add(node, m_ParentNode); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_ParentNode->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_ParentNode); } } if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if (m_ThreadIsRunning || !m_Visible) return; if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked())) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); DoFiberTracking(); } void QmitkStreamlineTrackingView::OnParameterChanged() { UpdateGui(); if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { if (m_FirstInteractiveRun) { QMessageBox::information(nullptr, "Information", "Place and move a spherical seed region anywhere in the image by left-clicking and dragging. If the seed region is colored red, tracking is in progress. If the seed region is colored white, tracking is finished.\nPlacing the seed region for the first time in a newly selected dataset might cause a short delay, since the tracker needs to be initialized."); m_FirstInteractiveRun = false; } QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); disconnect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } } void QmitkStreamlineTrackingView::Activated() { } void QmitkStreamlineTrackingView::Deactivated() { } void QmitkStreamlineTrackingView::Visible() { m_Visible = true; } void QmitkStreamlineTrackingView::Hidden() { m_Visible = false; m_Controls->m_InteractiveBox->setChecked(false); ToggleInteractive(); } void QmitkStreamlineTrackingView::OnSliceChanged() { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (!m_ThreadIsRunning && m_TrackingHandler != nullptr) { delete m_TrackingHandler; m_TrackingHandler = nullptr; m_DeleteTrackingHandler = false; + m_LastPrior = ""; + if (m_TrackingPriorHandler != nullptr) + delete m_TrackingPriorHandler; } else if (m_ThreadIsRunning) { m_DeleteTrackingHandler = true; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { std::vector< mitk::DataNode::Pointer > last_nodes = m_InputImageNodes; m_InputImageNodes.clear(); m_InputImages.clear(); m_AdditionalInputImages.clear(); bool retrack = false; for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData())) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr) { int dim = img->GetDimension(); unsigned int* dimensions = img->GetDimensions(); if (dim==4 && dimensions[3]%3==0) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if (dim==3) { m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } } } // sometimes the OnSelectionChanged event is sent twice and actually no selection has changed for the first event. We need to catch that. if (last_nodes.size() == m_InputImageNodes.size()) { bool same_nodes = true; for (unsigned int i=0; im_TensorImageLabel->setText("select in data-manager"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageBox->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestBox->setVisible(false); m_Controls->m_ForestLabel->setVisible(false); m_Controls->commandLinkButton->setEnabled(false); m_Controls->m_TrialsPerSeedBox->setEnabled(false); m_Controls->m_TrialsPerSeedLabel->setEnabled(false); m_Controls->m_TargetImageBox->setVisible(false); m_Controls->m_TargetImageLabel->setVisible(false); if (m_Controls->m_InteractiveBox->isChecked()) { m_Controls->m_InteractiveSeedingFrame->setVisible(true); m_Controls->m_StaticSeedingFrame->setVisible(false); m_Controls->commandLinkButton_2->setVisible(false); m_Controls->commandLinkButton->setVisible(false); } else { m_Controls->m_InteractiveSeedingFrame->setVisible(false); m_Controls->m_StaticSeedingFrame->setVisible(true); m_Controls->commandLinkButton_2->setVisible(m_ThreadIsRunning); m_Controls->commandLinkButton->setVisible(!m_ThreadIsRunning); } if (m_Controls->m_EpConstraintsBox->currentIndex()>0) { m_Controls->m_TargetImageBox->setVisible(true); m_Controls->m_TargetImageLabel->setVisible(true); } // trials per seed are only important for probabilistic tractography if (m_Controls->m_ModeBox->currentIndex()==1) { m_Controls->m_TrialsPerSeedBox->setEnabled(true); m_Controls->m_TrialsPerSeedLabel->setEnabled(true); } if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText( ( std::to_string(m_InputImageNodes.size()) + " images selected").c_str() ); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestBox->setVisible(true); m_Controls->m_ForestLabel->setVisible(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } } void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) { m_ThreadIsRunning = start; if (!m_Controls->m_InteractiveBox->isChecked()) { m_Controls->commandLinkButton_2->setVisible(start); m_Controls->commandLinkButton->setVisible(!start); m_Controls->m_InteractiveBox->setEnabled(!start); m_Controls->m_StatusTextBox->setVisible(start); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_ThreadIsRunning || m_InputImages.empty() || !m_Visible) return; if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) return; StartStopTrackingGui(true); m_Tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImages.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); StartStopTrackingGui(false); return; } -// if (m_FirstTensorProbRun) -// { -// QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. ODFs are sharpened (min-max normalized and raised to the power of 4). TEND parameters are ignored."); -// m_FirstTensorProbRun = false; -// } - if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::TensorImage::ItkTensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkImg ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); dynamic_cast(m_TrackingHandler)->SetIsOdfFromTensor(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (int i=0; i<(int)m_InputImages.size(); i++) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(i)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(m_TrackingHandler)->AddTensorImage(itkImg); } if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(m_TrackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); dynamic_cast(m_TrackingHandler)->SetOdfImage(itkImg); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestBox->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography (raw dMRI tractography) selected. Did you accidentally select the raw diffusion-weighted image in the datamanager?"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestBox->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { int num_previous_directions = (forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { int num_previous_directions = (forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); StartStopTrackingGui(false); return; } } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); StartStopTrackingGui(false); return; } try { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(0)); caster->SetCopyMemFlag(true); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(itkImg); } dynamic_cast(m_TrackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); } catch(...) { QMessageBox::information(nullptr, "Error", "Peak tracker could not be initialized. Is your input image in the correct format (4D float image, peaks in the 4th dimension)?"); StartStopTrackingGui(false); return; } } m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { m_Tracker->SetSeedPoints(m_SeedPoints); } else if (m_Controls->m_SeedImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetSeedImage(mask); } if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetMaskImage(mask); } if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetStoppingRegions(mask); } if (m_Controls->m_TargetImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetTargetRegions(mask); } if (m_Controls->m_PriorImageBox->GetSelectedNode().IsNotNull()) { - typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(dynamic_cast(m_Controls->m_PriorImageBox->GetSelectedNode()->GetData())); - caster->SetCopyMemFlag(true); - caster->Update(); - mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); - mitk::TrackingDataHandler* trackingPriorHandler = new mitk::TrackingHandlerPeaks(); - dynamic_cast(trackingPriorHandler)->SetPeakImage(itkImg); - dynamic_cast(trackingPriorHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); - - trackingPriorHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); - trackingPriorHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); - trackingPriorHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); - trackingPriorHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); - trackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); - - m_Tracker->SetTrackingPriorHandler(trackingPriorHandler); + if (m_LastPrior!=m_Controls->m_PriorImageBox->GetSelectedNode()->GetUID() || m_TrackingPriorHandler==nullptr) + { + typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(dynamic_cast(m_Controls->m_PriorImageBox->GetSelectedNode()->GetData())); + caster->SetCopyMemFlag(true); + caster->Update(); + mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); + m_TrackingPriorHandler = new mitk::TrackingHandlerPeaks(); + dynamic_cast(m_TrackingPriorHandler)->SetPeakImage(itkImg); + dynamic_cast(m_TrackingPriorHandler)->SetPeakThreshold(0.0); + m_LastPrior = m_Controls->m_PriorImageBox->GetSelectedNode()->GetUID(); + } + + m_TrackingPriorHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); + m_TrackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); + + m_Tracker->SetTrackingPriorHandler(m_TrackingPriorHandler); m_Tracker->SetTrackingPriorWeight(m_Controls->m_PriorWeightBox->value()); m_Tracker->SetTrackingPriorAsMask(m_Controls->m_PriorAsMaskBox->isChecked()); m_Tracker->SetIntroduceDirectionsFromPrior(m_Controls->m_NewDirectionsFromPriorBox->isChecked()); } + else if (m_Controls->m_PriorImageBox->GetSelectedNode().IsNull()) + m_Tracker->SetTrackingPriorHandler(nullptr); if (m_Controls->m_ExclusionImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_ExclusionImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetExclusionRegions(mask); } // Endpoint constraints switch (m_Controls->m_EpConstraintsBox->currentIndex()) { case 0: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); m_Tracker->SetTargetRegions(nullptr); break; case 1: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); break; case 2: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); break; case 3: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); break; case 4: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); break; case 5: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); break; case 6: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); break; } if (m_Tracker->GetEndpointConstraint()!=itk::StreamlineTrackingFilter::EndpointConstraints::NONE && m_Controls->m_TargetImageBox->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Error", "Endpoint constraints are used but no target image is set!"); StartStopTrackingGui(false); return; } else if (m_Tracker->GetEndpointConstraint()==itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET && (m_Controls->m_SeedImageBox->GetSelectedNode().IsNull()|| m_Controls->m_TargetImageBox->GetSelectedNode().IsNull()) ) { QMessageBox::information(nullptr, "Error", "Endpoint constraint EPS_IN_SEED_AND_TARGET is used but no target or no seed image is set!"); StartStopTrackingGui(false); return; } m_Tracker->SetInterpolateMasks(m_Controls->m_MaskInterpolationBox->isChecked()); m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); m_Tracker->SetTrialsPerSeed(m_Controls->m_TrialsPerSeedBox->value()); m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); m_Tracker->SetLoopCheck(m_Controls->m_LoopCheckBox->value()); m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); m_ParentNode = m_InputImageNodes.at(0); m_TrackingThread.start(QThread::LowestPriority); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h index 0f6db98e80..d5d0d6ee99 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h @@ -1,151 +1,153 @@ /*=================================================================== 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 QmitkStreamlineTrackingView_h #define QmitkStreamlineTrackingView_h #include #include "ui_QmitkStreamlineTrackingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include class QmitkStreamlineTrackingView; class QmitkStreamlineTrackingWorker : public QObject { Q_OBJECT public: QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view); public slots: void run(); private: QmitkStreamlineTrackingView* m_View; }; /*! \brief View for tensor based deterministic streamline fiber tracking. */ class QmitkStreamlineTrackingView : public QmitkAbstractView, public mitk::ILifecycleAwarePart { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; typedef itk::Image< unsigned int, 3 > ItkUintImgType; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; typedef itk::Image< float, 3 > ItkFloatImageType; typedef itk::StreamlineTrackingFilter TrackerType; QmitkStreamlineTrackingView(); virtual ~QmitkStreamlineTrackingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; TrackerType::Pointer m_Tracker; QmitkStreamlineTrackingWorker m_TrackingWorker; QThread m_TrackingThread; virtual void Activated() override; virtual void Deactivated() override; virtual void Visible() override; virtual void Hidden() override; protected slots: void DoFiberTracking(); ///< start fiber tracking void UpdateGui(); void ToggleInteractive(); void DeleteTrackingHandler(); void OnParameterChanged(); void InteractiveSeedChanged(bool posChanged=false); void ForestSwitched(); void OutputStyleSwitched(); void AfterThread(); ///< update gui etc. after tracking has finished void BeforeThread(); ///< start timer etc. void TimerUpdate(); void StopTractography(); void OnSliceChanged(); protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkStreamlineTrackingViewControls* m_Controls; protected slots: private: void StartStopTrackingGui(bool start); std::vector< itk::Point > m_SeedPoints; mitk::DataNode::Pointer m_ParentNode; mitk::DataNode::Pointer m_InteractiveNode; mitk::DataNode::Pointer m_InteractivePointSetNode; std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input image nodes std::vector< mitk::Image::ConstPointer > m_InputImages; ///< input images std::vector< mitk::Image::ConstPointer > m_AdditionalInputImages; bool m_FirstTensorProbRun; bool m_FirstInteractiveRun; mitk::TrackingDataHandler* m_TrackingHandler; bool m_ThreadIsRunning; QTimer* m_TrackingTimer; bool m_DeleteTrackingHandler; QmitkSliceNavigationListener m_SliceChangeListener; bool m_Visible; + mitk::Identifiable::UIDType m_LastPrior; + mitk::TrackingDataHandler* m_TrackingPriorHandler; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED