diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp index ea2eac7eb8..50d790be8f 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp @@ -1,357 +1,296 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include template int StartPeakExtraction(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("image", "i", mitkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false); parser.addArgument("outroot", "o", mitkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false); - parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask", "mask image"); - parser.addArgument("normalization", "n", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true); - parser.addArgument("numpeaks", "p", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true); - parser.addArgument("peakthres", "r", mitkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true); - parser.addArgument("abspeakthres", "a", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true); - parser.addArgument("shConvention", "s", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", std::string("MITK"), true); - parser.addArgument("noFlip", "f", mitkCommandLineParser::Bool, "No flip", "do not flip input image to match MITK coordinate convention"); - parser.addArgument("clusterThres", "c", mitkCommandLineParser::Float, "Clustering threshold", "directions closer together than the specified angular threshold will be clustered (in rad)", 0.9); - parser.addArgument("flipX", "fx", mitkCommandLineParser::Bool, "Flip X", "Flip peaks in x direction"); - parser.addArgument("flipY", "fy", mitkCommandLineParser::Bool, "Flip Y", "Flip peaks in y direction"); - parser.addArgument("flipZ", "fz", mitkCommandLineParser::Bool, "Flip Z", "Flip peaks in z direction"); + parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "Mask", "mask image"); + parser.addArgument("normalization", "", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true); + parser.addArgument("numpeaks", "", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true); + + parser.addArgument("rel_peakthr", "", mitkCommandLineParser::Float, "Relative peak threshold", "peak threshold relative to largest peak", 0.4, true); + parser.addArgument("abs_peakthr", "", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak magnitude threshold", 0.03, true); + parser.addArgument("angular_thr", "", mitkCommandLineParser::Float, "Angular threshold", "in degree", 15); + + parser.addArgument("shConvention", "", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MRtrix, FSL)", std::string("MRtrix"), true); + + + parser.addArgument("flipX", "", mitkCommandLineParser::Bool, "Flip X", "Flip peaks in x direction"); + parser.addArgument("flipY", "", mitkCommandLineParser::Bool, "Flip Y", "Flip peaks in y direction"); + parser.addArgument("flipZ", "", mitkCommandLineParser::Bool, "Flip Z", "Flip peaks in z direction"); + parser.addArgument("scale_by_gfa", "", mitkCommandLineParser::Bool, "Scale by GFA", "Scale ODF values and peaks by GFA"); + parser.setCategory("Preprocessing Tools"); parser.setTitle("Peak Extraction"); parser.setDescription(""); parser.setContributor("MIC"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["image"]); std::string outRoot = us::any_cast(parsedArgs["outroot"]); // optional arguments std::string maskImageName(""); if (parsedArgs.count("mask")) maskImageName = us::any_cast(parsedArgs["mask"]); int normalization = 1; if (parsedArgs.count("normalization")) normalization = us::any_cast(parsedArgs["normalization"]); int numPeaks = 2; if (parsedArgs.count("numpeaks")) numPeaks = us::any_cast(parsedArgs["numpeaks"]); - float peakThres = 0.4; - if (parsedArgs.count("peakthres")) - peakThres = us::any_cast(parsedArgs["peakthres"]); + float rel_peakthr = 0.4; + if (parsedArgs.count("rel_peakthr")) + rel_peakthr = us::any_cast(parsedArgs["rel_peakthr"]); - float absPeakThres = 0.06; - if (parsedArgs.count("abspeakthres")) - absPeakThres = us::any_cast(parsedArgs["abspeakthres"]); + float abs_peakthr = 0.03; + if (parsedArgs.count("abs_peakthr")) + abs_peakthr = us::any_cast(parsedArgs["abs_peakthr"]); - float clusterThres = 0.9; - if (parsedArgs.count("clusterThres")) - clusterThres = us::any_cast(parsedArgs["clusterThres"]); + float angular_thr = 15; + if (parsedArgs.count("angular_thr")) + angular_thr = us::any_cast(parsedArgs["angular_thr"]); + angular_thr = cos((float)angular_thr*itk::Math::pi / 180); - bool noFlip = false; - if (parsedArgs.count("noFlip")) - noFlip = us::any_cast(parsedArgs["noFlip"]); + bool scale_by_gfa = false; + if (parsedArgs.count("scale_by_gfa")) + scale_by_gfa = us::any_cast(parsedArgs["scale_by_gfa"]); bool flipX = false; if (parsedArgs.count("flipX")) flipX = us::any_cast(parsedArgs["flipX"]); bool flipY = false; if (parsedArgs.count("flipY")) flipY = us::any_cast(parsedArgs["flipY"]); bool flipZ = false; if (parsedArgs.count("flipZ")) flipZ = us::any_cast(parsedArgs["flipZ"]); std::cout << "image: " << imageName; std::cout << "outroot: " << outRoot; if (!maskImageName.empty()) std::cout << "mask: " << maskImageName; else std::cout << "no mask image selected"; std::cout << "numpeaks: " << numPeaks; - std::cout << "peakthres: " << peakThres; - std::cout << "abspeakthres: " << absPeakThres; + std::cout << "peakthres: " << rel_peakthr; + std::cout << "abspeakthres: " << abs_peakthr; std::cout << "shOrder: " << shOrder; try { mitk::Image::Pointer image = mitk::IOUtil::Load(imageName); mitk::Image::Pointer mask = mitk::IOUtil::Load(maskImageName); typedef itk::Image ItkUcharImgType; - typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; + typedef itk::OdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; typename MaximaExtractionFilterType::Pointer peak_extraction_filter = MaximaExtractionFilterType::New(); - int toolkitConvention = 0; - - if (parsedArgs.count("shConvention")) - { - std::string convention = us::any_cast(parsedArgs["shConvention"]).c_str(); - if ( boost::algorithm::equals(convention, "FSL") ) - { - toolkitConvention = 1; - std::cout << "Using FSL SH-basis"; - } - else if ( boost::algorithm::equals(convention, "MRtrix") ) - { - toolkitConvention = 2; - std::cout << "Using MRtrix SH-basis"; - } - else - std::cout << "Using MITK SH-basis"; - } - else - std::cout << "Using MITK SH-basis"; ItkUcharImgType::Pointer itkMaskImage = nullptr; if (mask.IsNotNull()) { try{ itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(mask, itkMaskImage); peak_extraction_filter->SetMaskImage(itkMaskImage); } catch(...) { } } - if (toolkitConvention>0) + if (parsedArgs.count("shConvention")) { - std::cout << "Converting coefficient image to MITK format"; - typedef itk::ShCoefficientImageImporter< float, shOrder > ConverterType; - typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(image); - caster->Update(); - itk::Image< float, 4 >::Pointer itkImage = caster->GetOutput(); - - typename ConverterType::Pointer converter = ConverterType::New(); - - if (noFlip) - { - converter->SetInputImage(itkImage); - } - else - { - std::cout << "Flipping image"; - itk::FixedArray flipAxes; - flipAxes[0] = true; - flipAxes[1] = true; - flipAxes[2] = false; - flipAxes[3] = false; - itk::FlipImageFilter< itk::Image< float, 4 > >::Pointer flipper = itk::FlipImageFilter< itk::Image< float, 4 > >::New(); - flipper->SetInput(itkImage); - flipper->SetFlipAxes(flipAxes); - flipper->Update(); - itk::Image< float, 4 >::Pointer flipped = flipper->GetOutput(); - itk::Matrix< double,4,4 > m = itkImage->GetDirection(); m[0][0] *= -1; m[1][1] *= -1; - flipped->SetDirection(m); - - itk::Point< float, 4 > o = itkImage->GetOrigin(); - o[0] -= (flipped->GetLargestPossibleRegion().GetSize(0)-1); - o[1] -= (flipped->GetLargestPossibleRegion().GetSize(1)-1); - flipped->SetOrigin(o); - converter->SetInputImage(flipped); - } - - std::cout << "Starting conversion"; - switch (toolkitConvention) - { - case 1: + std::string convention = us::any_cast(parsedArgs["shConvention"]).c_str(); + if ( convention=="FSL" ) peak_extraction_filter->SetToolkit(MaximaExtractionFilterType::FSL); - break; - case 2: + else peak_extraction_filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); - break; - default: - peak_extraction_filter->SetToolkit(MaximaExtractionFilterType::FSL); - break; - } - converter->GenerateData(); - peak_extraction_filter->SetInput(converter->GetCoefficientImage()); } else + peak_extraction_filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); + + try{ + typedef mitk::ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; + typename CasterType::Pointer caster = CasterType::New(); + caster->SetInput(image); + caster->Update(); + peak_extraction_filter->SetInput(caster->GetOutput()); + } + catch(...) { - try{ - typedef mitk::ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; - typename CasterType::Pointer caster = CasterType::New(); - caster->SetInput(image); - caster->Update(); - peak_extraction_filter->SetInput(caster->GetOutput()); - } - catch(...) - { - std::cout << "wrong image type"; - return EXIT_FAILURE; - } + std::cout << "wrong image type"; + return EXIT_FAILURE; } peak_extraction_filter->SetMaxNumPeaks(numPeaks); - peak_extraction_filter->SetPeakThreshold(peakThres); - peak_extraction_filter->SetAbsolutePeakThreshold(absPeakThres); - peak_extraction_filter->SetAngularThreshold(1); - peak_extraction_filter->SetClusteringThreshold(clusterThres); + peak_extraction_filter->SetRelativePeakThreshold(rel_peakthr); + peak_extraction_filter->SetAbsolutePeakThreshold(abs_peakthr); + peak_extraction_filter->SetAngularThreshold(angular_thr); peak_extraction_filter->SetFlipX(flipX); peak_extraction_filter->SetFlipY(flipY); peak_extraction_filter->SetFlipZ(flipZ); + peak_extraction_filter->SetScaleByGfa(scale_by_gfa); switch (normalization) { case 0: peak_extraction_filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: peak_extraction_filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: peak_extraction_filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } std::cout << "Starting extraction"; peak_extraction_filter->Update(); // write direction image { typename MaximaExtractionFilterType::PeakImageType::Pointer itkImg = peak_extraction_filter->GetPeakImage(); std::string outfilename = outRoot; outfilename.append("_PEAKS.nrrd"); typedef itk::ImageFileWriter< typename MaximaExtractionFilterType::PeakImageType > WriterType; typename WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(itkImg); writer->Update(); } // write num directions image { ItkUcharImgType::Pointer numDirImage = peak_extraction_filter->GetNumDirectionsImage(); if (itkMaskImage.IsNotNull()) { numDirImage->SetDirection(itkMaskImage->GetDirection()); numDirImage->SetOrigin(itkMaskImage->GetOrigin()); } std::string outfilename = outRoot.c_str(); outfilename.append("_NUM_PEAKS.nrrd"); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(numDirImage); writer->Update(); } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } /*! \brief Extract maxima in the input spherical harmonics image. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("image", "i", mitkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false); parser.addArgument("shOrder", "sh", mitkCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order"); parser.addArgument("outroot", "o", mitkCommandLineParser::OutputDirectory, "Output directory", "output root", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask", "mask image"); parser.addArgument("normalization", "n", mitkCommandLineParser::Int, "Normalization", "0=no norm, 1=max norm, 2=single vec norm", 1, true); parser.addArgument("numpeaks", "p", mitkCommandLineParser::Int, "Max. number of peaks", "maximum number of extracted peaks", 2, true); parser.addArgument("peakthres", "r", mitkCommandLineParser::Float, "Peak threshold", "peak threshold relative to largest peak", 0.4, true); parser.addArgument("abspeakthres", "a", mitkCommandLineParser::Float, "Absolute peak threshold", "absolute peak threshold weighted with local GFA value", 0.06, true); parser.addArgument("shConvention", "s", mitkCommandLineParser::String, "Use specified SH-basis", "use specified SH-basis (MITK, FSL, MRtrix)", std::string("MITK"), true); parser.addArgument("noFlip", "f", mitkCommandLineParser::Bool, "No flip", "do not flip input image to match MITK coordinate convention"); parser.setCategory("Preprocessing Tools"); parser.setTitle("Peak Extraction"); parser.setDescription("Extract maxima in the input spherical harmonics image."); parser.setContributor("MIC"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; int shOrder = -1; if (parsedArgs.count("shOrder")) shOrder = us::any_cast(parsedArgs["shOrder"]); switch (shOrder) { case 4: return StartPeakExtraction<4>(argc, argv); case 6: return StartPeakExtraction<6>(argc, argv); case 8: return StartPeakExtraction<8>(argc, argv); case 10: return StartPeakExtraction<10>(argc, argv); case 12: return StartPeakExtraction<12>(argc, argv); } return EXIT_FAILURE; } diff --git a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp index ebf019b60f..35f94b63e2 100644 --- a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkDiffusionCoreIOMimeTypes.cpp @@ -1,642 +1,646 @@ /*=================================================================== 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 itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); // Simple NRRD files should only be considered for this mime type if they contain // corresponding tags if( io->CanReadFile(path.c_str())) { 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 { + std::string ext = itksys::SystemTools::GetFilenameExtension(path); + if (ext==".shi") + 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) { 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/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake index 5a2471e704..b8972e792d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake @@ -1,163 +1,162 @@ set(CPP_FILES # DicomImport # DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp DicomImport/mitkDiffusionDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp DicomImport/mitkDiffusionHeaderSiemensMosaicDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.cpp # Properties IODataStructures/Properties/mitkBValueMapProperty.cpp IODataStructures/Properties/mitkGradientDirectionsProperty.cpp IODataStructures/Properties/mitkMeasurementFrameProperty.cpp IODataStructures/Properties/mitkDiffusionPropertyHelper.cpp IODataStructures/Properties/mitkNodePredicateIsDWI.cpp # Serializer IODataStructures/Properties/mitkBValueMapPropertySerializer.cpp IODataStructures/Properties/mitkGradientDirectionsPropertySerializer.cpp IODataStructures/Properties/mitkMeasurementFramePropertySerializer.cpp # DataStructures -> Odf IODataStructures/OdfImages/mitkOdfImageSource.cpp IODataStructures/OdfImages/mitkOdfImage.cpp IODataStructures/mitkShImage.cpp IODataStructures/mitkShImageSource.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImage.cpp # DataStructures -> Peaks IODataStructures/mitkPeakImage.cpp Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/itkDwiGradientLengthCorrectionFilter.cpp # Registration Algorithms & Co. Algorithms/Registration/mitkRegistrationWrapper.cpp Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp # Algorithms/Registration/mitkRegistrationMethodITK4.cpp Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp # MultishellProcessing Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection include/mitkDiffusionFunctionCollection.h # Rendering include/Rendering/mitkOdfVtkMapper2D.h # Reconstruction include/Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h include/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h include/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h include/Algorithms/Reconstruction/itkPointShell.h include/Algorithms/Reconstruction/itkOrientationDistributionFunction.h include/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.h include/Algorithms/Reconstruction/itkMultiTensorImageFilter.h # Fitting functions include/Algorithms/Reconstruction/FittingFunctions/mitkAbstractFitter.h include/Algorithms/Reconstruction/FittingFunctions/mitkMultiTensorFitter.h include/Algorithms/Reconstruction/FittingFunctions/mitkBallStickFitter.h # MultishellProcessing include/Algorithms/Reconstruction/MultishellProcessing/itkRadialMultishellToSingleshellImageFilter.h include/Algorithms/Reconstruction/MultishellProcessing/itkDWIVoxelFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.h # Properties include/IODataStructures/Properties/mitkBValueMapProperty.h include/IODataStructures/Properties/mitkGradientDirectionsProperty.h include/IODataStructures/Properties/mitkMeasurementFrameProperty.h include/IODataStructures/Properties/mitkDiffusionPropertyHelper.h include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.h # Algorithms include/Algorithms/itkDiffusionOdfGeneralizedFaImageFilter.h include/Algorithms/itkDiffusionOdfPrepareVisualizationImageFilter.h include/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h include/Algorithms/itkTensorDerivedMeasurementsFilter.h include/Algorithms/itkBrainMaskExtractionImageFilter.h include/Algorithms/itkB0ImageExtractionImageFilter.h include/Algorithms/itkB0ImageExtractionToSeparateImageFilter.h include/Algorithms/itkTensorImageToDiffusionImageFilter.h include/Algorithms/itkTensorToL2NormImageFilter.h include/Algorithms/itkGaussianInterpolateImageFunction.h include/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h include/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h include/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h include/Algorithms/itkCartesianToPolarVectorImageFilter.h include/Algorithms/itkPolarToCartesianVectorImageFilter.h include/Algorithms/itkDistanceMapFilter.h include/Algorithms/itkProjectionFilter.h include/Algorithms/itkResidualImageFilter.h include/Algorithms/itkExtractChannelFromRgbaImageFilter.h include/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h include/Algorithms/itkMergeDiffusionImagesFilter.h - include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h include/Algorithms/itkShCoefficientImageImporter.h include/Algorithms/itkShCoefficientImageExporter.h include/Algorithms/itkOdfMaximaExtractionFilter.h include/Algorithms/itkResampleDwiImageFilter.h include/Algorithms/itkDwiGradientLengthCorrectionFilter.h include/Algorithms/itkAdcImageFilter.h include/Algorithms/itkDwiNormilzationFilter.h include/Algorithms/itkSplitDWImageFilter.h include/Algorithms/itkRemoveDwiChannelFilter.h include/Algorithms/itkExtractDwiChannelFilter.h include/Algorithms/itkFlipPeaksFilter.h include/Algorithms/itkShToOdfImageFilter.h include/Algorithms/itkFourDToVectorImageFilter.h include/Algorithms/itkVectorImageToFourDImageFilter.h include/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h include/Algorithms/itkNonLocalMeansDenoisingFilter.h include/Algorithms/itkVectorImageToImageFilter.h include/Algorithms/itkSplitVectorImageFilter.h ) set( TOOL_FILES ) diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp deleted file mode 100644 index 7bec5289de..0000000000 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp +++ /dev/null @@ -1,430 +0,0 @@ -/*=================================================================== - -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 __itkFiniteDiffOdfMaximaExtractionFilter_cpp -#define __itkFiniteDiffOdfMaximaExtractionFilter_cpp - -#include "itkFiniteDiffOdfMaximaExtractionFilter.h" -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include - -using namespace boost::math; - -namespace itk { - - -static bool CompareVectors(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) -{ - return (v1.magnitude()>v2.magnitude()); -} - -template< class PixelType, int ShOrder, int NrOdfDirections > -FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> -::FiniteDiffOdfMaximaExtractionFilter() - : m_NormalizationMethod(MAX_VEC_NORM) - , m_MaxNumPeaks(2) - , m_PeakThreshold(0.4) - , m_AbsolutePeakThreshold(0) - , m_ClusteringThreshold(0.9) - , m_AngularThreshold(0.7) - , m_NumCoeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder) - , m_Toolkit(MRTRIX) - , m_ApplyDirectionMatrix(false) -{ - this->SetNumberOfRequiredInputs(1); -} - -template< class PixelType, int ShOrder, int NrOdfDirections > -void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> -::FindCandidatePeaks(OdfType& odf, double thr, std::vector< DirectionType >& container) -{ - double gfa = odf.GetGeneralizedFractionalAnisotropy(); - //Find the peaks using a finite difference method - bool flag = true; - vnl_vector_fixed< bool, NrOdfDirections > used; used.fill(false); - //Find the peaks - for (int i=0; ithr && val*gfa>m_AbsolutePeakThreshold) // limit to one hemisphere ??? - { - flag = true; - std::vector< int > neighbours = odf.GetNeighbors(i); - for (unsigned int j=0; j -std::vector< vnl_vector_fixed< double, 3 > > FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>::MeanShiftClustering(std::vector< DirectionType >& inDirs) -{ - std::vector< DirectionType > outDirs; - if (inDirs.empty()) - return inDirs; - - DirectionType oldMean, currentMean, workingMean; - std::vector< int > touched; - - // initialize - touched.resize(inDirs.size(), 0); - bool free = true; - currentMean = inDirs[0]; // initialize first seed - while (free) - { - oldMean.fill(0.0); - - // start mean-shift clustering - float angle = 0.0; - int counter = 0; - while ((currentMean-oldMean).magnitude()>0.0001) - { - counter = 0; - oldMean = currentMean; - workingMean = oldMean; - workingMean.normalize(); - currentMean.fill(0.0); - for (unsigned int i=0; i=m_ClusteringThreshold) - { - currentMean += inDirs[i]; - touched[i] = 1; - counter++; - } - else if (-angle>=m_ClusteringThreshold) - { - currentMean -= inDirs[i]; - touched[i] = 1; - counter++; - } - } - } - - // found stable mean - if (counter>0) - { - float mag = currentMean.magnitude(); - if (mag>0) - { - currentMean /= mag; - outDirs.push_back(currentMean); - } - } - - // find next unused seed - free = false; - for (unsigned int i=0; i -void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> -::BeforeThreadedGenerateData() -{ - typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); - itk::Vector spacing = ShCoeffImage->GetSpacing(); - double minSpacing = spacing[0]; - if (spacing[1]GetOrigin(); - itk::Matrix direction = ShCoeffImage->GetDirection(); - ImageRegion<3> imageRegion = ShCoeffImage->GetLargestPossibleRegion(); - - if (m_MaskImage.IsNotNull()) - { - origin = m_MaskImage->GetOrigin(); - direction = m_MaskImage->GetDirection(); - imageRegion = m_MaskImage->GetLargestPossibleRegion(); - } - - itk::Vector spacing3 = ShCoeffImage->GetSpacing(); - itk::Point origin3 = ShCoeffImage->GetOrigin(); - itk::Matrix direction3 = ShCoeffImage->GetDirection(); - itk::ImageRegion<3> imageRegion3 = ShCoeffImage->GetLargestPossibleRegion(); - - itk::Vector spacing4; - itk::Point origin4; - itk::Matrix direction4; - itk::ImageRegion<4> imageRegion4; - - spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; - origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin4[3] = 0; - for (int r=0; r<3; r++) - for (int c=0; c<3; c++) - direction4[r][c] = direction3[r][c]; - direction4[3][3] = 1; - imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); - imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); - imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); - imageRegion4.SetSize(3, m_MaxNumPeaks*3); - - m_PeakImage = PeakImageType::New(); - m_PeakImage->SetSpacing( spacing4 ); - m_PeakImage->SetOrigin( origin4 ); - m_PeakImage->SetDirection( direction4 ); - m_PeakImage->SetRegions( imageRegion4 ); - m_PeakImage->Allocate(); - m_PeakImage->FillBuffer(0.0); - - if (m_MaskImage.IsNull()) - { - m_MaskImage = ItkUcharImgType::New(); - m_MaskImage->SetSpacing( spacing ); - m_MaskImage->SetOrigin( origin ); - m_MaskImage->SetDirection( direction ); - m_MaskImage->SetRegions( imageRegion ); - m_MaskImage->Allocate(); - m_MaskImage->FillBuffer(1); - } - m_NumDirectionsImage = ItkUcharImgType::New(); - m_NumDirectionsImage->SetSpacing( spacing ); - m_NumDirectionsImage->SetOrigin( origin ); - m_NumDirectionsImage->SetDirection( direction ); - m_NumDirectionsImage->SetRegions( imageRegion ); - m_NumDirectionsImage->Allocate(); - m_NumDirectionsImage->FillBuffer(0); - - // calculate SH basis - OdfType odf; - vnl_matrix< double > sphCoords; - std::vector< DirectionType > dirs; - for (int i=0; i -void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> -::AfterThreadedGenerateData() -{ - -} - -template< class PixelType, int ShOrder, int NrOdfDirections > -void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> -::ThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, ThreadIdType threadID ) -{ - typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); - - ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, outputRegionForThread ); - - OdfType odf; - while( !cit.IsAtEnd() ) - { - typename CoefficientImageType::IndexType idx3 = cit.GetIndex(); - if (m_MaskImage->GetPixel(idx3)==0) - { - ++cit; - continue; - } - - CoefficientPixelType c = cit.Get(); - - // calculate ODF - double max = 0; - odf.Fill(0.0); - for (int i=0; imax) - max = odf[i]; - } - if (max<0.0001) - { - ++cit; - continue; - } - - std::vector< DirectionType > candidates, peaks, temp; - peaks.clear(); - max *= m_PeakThreshold; // relative threshold - FindCandidatePeaks(odf, max, candidates); // find all local maxima - candidates = MeanShiftClustering(candidates); // cluster maxima - - vnl_matrix sphCoords; - CreateDirMatrix(candidates, sphCoords); // convert candidate peaks to spherical angles - vnl_matrix< float > shBasis; - if (m_Toolkit==Toolkit::MRTRIX) - shBasis = mitk::sh::CalcShBasisForDirections(ShOrder, sphCoords); - else - shBasis = mitk::sh::CalcShBasisForDirections(ShOrder, sphCoords, false); - - max = 0.0; - for (unsigned int i=0; imax) - max = val; - peaks.push_back(candidates[i]*val); - } - std::sort( peaks.begin(), peaks.end(), CompareVectors ); // sort peaks - - // kick out directions to close to a larger direction (too far away to cluster but too close to keep) - unsigned int m = peaks.size(); - if ( m>m_MaxNumPeaks ) - m = m_MaxNumPeaks; - for (unsigned int i=0; im_AngularThreshold && val idx4; idx4[0] = idx3[0]; idx4[1] = idx3[1]; idx4[2] = idx3[2]; - - // fill output image - unsigned int num = peaks.size(); - if ( num>m_MaxNumPeaks ) - num = m_MaxNumPeaks; - for (unsigned int i=0; iGetDirection()*dir; - - if (m_FlipX) - dir[0] = -dir[0]; - if (m_FlipY) - dir[1] = -dir[1]; - if (m_FlipZ) - dir[2] = -dir[2]; - - for (unsigned int j = 0; j<3; j++) - { - idx4[3] = i*3 + j; - m_PeakImage->SetPixel(idx4, dir[j]); - } - } - m_NumDirectionsImage->SetPixel(idx3, num); - ++cit; - } - MITK_INFO << "Thread " << threadID << " finished extraction"; -} - -// convert cartesian to spherical coordinates -template< class PixelType, int ShOrder, int NrOdfDirections > -void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> -::CreateDirMatrix(const std::vector< DirectionType >& dir, vnl_matrix& sphCoords) -{ - sphCoords.set_size(3, dir.size()); - for (unsigned int i=0; i -#include - -namespace itk{ - -/** -* \brief Extract ODF peaks by searching all local maxima on a densely sampled ODF und clustering these maxima to get the underlying fiber direction. -* NrOdfDirections: number of sampling points on the ODF surface (about 20000 is a good value) -*/ - -template< class PixelType, int ShOrder, int NrOdfDirections > -class FiniteDiffOdfMaximaExtractionFilter : public ImageToImageFilter< Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, -Image< unsigned char, 3 > > -{ - - public: - - enum Toolkit { ///< SH coefficient convention (depends on toolkit) - FSL, - MRTRIX - }; - - enum NormalizationMethods { - NO_NORM, ///< no length normalization of the output peaks - SINGLE_VEC_NORM, ///< normalize the single peaks to length 1 - MAX_VEC_NORM ///< normalize all peaks according to their length in comparison to the largest peak in the respective voxel (0-1) - }; - - typedef FiniteDiffOdfMaximaExtractionFilter Self; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - typedef ImageToImageFilter< Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, - Image< unsigned char, 3 > > Superclass; - - /** Method for creation through the object factory. */ - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - - /** Runtime information support. */ - itkTypeMacro(FiniteDiffOdfMaximaExtractionFilter, ImageToImageFilter) - - typedef typename Superclass::InputImageType CoefficientImageType; - typedef typename CoefficientImageType::RegionType CoefficientImageRegionType; - typedef typename CoefficientImageType::PixelType CoefficientPixelType; - typedef typename Superclass::OutputImageType OutputImageType; - typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - typedef typename Superclass::InputImageRegionType InputImageRegionType; - typedef Image< float, 4 > PeakImageType; - - typedef OrientationDistributionFunction OdfType; - typedef itk::Image ItkUcharImgType; - - typedef vnl_vector_fixed< double, 3 > DirectionType; - - // input - itkSetMacro( MaxNumPeaks, unsigned int) ///< maximum number of peaks per voxel. if more peaks are detected, only the largest are kept. - itkSetMacro( PeakThreshold, double) ///< threshold on the peak length relative to the largest peak inside the current voxel - itkSetMacro( AbsolutePeakThreshold, double) ///< hard threshold on the peak length of all local maxima - itkSetMacro( ClusteringThreshold, double) ///< directions closer together than the specified angular threshold will be clustered (in rad) - itkSetMacro( AngularThreshold, double) ///< directions closer together than the specified threshold that remain after clustering are discarded (largest is kept) (in rad) - itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< only voxels inside the binary mask are processed - itkSetMacro( NormalizationMethod, NormalizationMethods) ///< normalization method of ODF peaks - itkSetMacro( FlipX, bool) ///< flip peaks in x direction - itkSetMacro( FlipY, bool) ///< flip peaks in y direction - itkSetMacro( FlipZ, bool) ///< flip peaks in z direction - itkSetMacro( ApplyDirectionMatrix, bool) - - // output - itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer ) - itkGetMacro( PeakImage, PeakImageType::Pointer ) - - itkSetMacro( Toolkit, Toolkit) ///< define SH coefficient convention (depends on toolkit) - itkGetMacro( Toolkit, Toolkit) ///< SH coefficient convention (depends on toolkit) - - protected: - FiniteDiffOdfMaximaExtractionFilter(); - ~FiniteDiffOdfMaximaExtractionFilter(){} - - void BeforeThreadedGenerateData(); - void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadID ); - void AfterThreadedGenerateData(); - - /** Extract all local maxima from the densely sampled ODF surface. Thresholding possible. **/ - void FindCandidatePeaks(OdfType& odf, double odfMax, std::vector< DirectionType >& inDirs); - - /** Cluster input directions within a certain angular threshold **/ - std::vector< DirectionType > MeanShiftClustering(std::vector< DirectionType >& inDirs); - - /** Convert direction vector to matrix **/ - void CreateDirMatrix(const std::vector< DirectionType >& dir, vnl_matrix& sphCoords); - - private: - - NormalizationMethods m_NormalizationMethod; ///< normalization method of ODF peaks - unsigned int m_MaxNumPeaks; ///< maximum number of peaks per voxel. if more peaks are detected, only the largest are kept. - double m_PeakThreshold; ///< threshold on the peak length relative to the largest peak inside the current voxel - double m_AbsolutePeakThreshold;///< hard threshold on the peak length of all local maxima - vnl_matrix< float > m_ShBasis; ///< container for evaluated SH base functions - double m_ClusteringThreshold; ///< directions closer together than the specified angular threshold will be clustered (in rad) - double m_AngularThreshold; ///< directions closer together than the specified threshold that remain after clustering are discarded (largest is kept) (in rad) - const int m_NumCoeffs; ///< number of spherical harmonics coefficients - - PeakImageType::Pointer m_PeakImage; - ItkUcharImgType::Pointer m_NumDirectionsImage; ///< number of peaks per voxel - ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed - - Toolkit m_Toolkit; - bool m_FlipX; - bool m_FlipY; - bool m_FlipZ; - bool m_ApplyDirectionMatrix; -}; - -} - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkFiniteDiffOdfMaximaExtractionFilter.cpp" -#endif - -#endif //__itkFiniteDiffOdfMaximaExtractionFilter_h_ - diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.cpp index 0a8bc8f26a..7a1f49a791 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.cpp @@ -1,636 +1,358 @@ /*=================================================================== 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 __itkOdfMaximaExtractionFilter_cpp #define __itkOdfMaximaExtractionFilter_cpp - #include "itkOdfMaximaExtractionFilter.h" +#include +#include #include +#include +#include #include #include #include #include #include #include -#include #include -#include -#include +#include using namespace boost::math; namespace itk { -template< class TOdfPixelType > -OdfMaximaExtractionFilter< TOdfPixelType >::OdfMaximaExtractionFilter(): - m_NormalizationMethod(MAX_VEC_NORM), - m_PeakThreshold(0.2), - m_MaxNumPeaks(10), - m_ShCoeffImage(nullptr), - m_OutputFiberBundle(nullptr), - m_NumDirectionsImage(nullptr), - m_DirectionImageContainer(nullptr) -{ +static bool CompareVectors(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) +{ + return (v1.magnitude()>v2.magnitude()); } -// solve ax? + bx? + cx + d = 0 using cardanos method -template< class TOdfPixelType > -bool OdfMaximaExtractionFilter::ReconstructQballImage() +template< class PixelType, int ShOrder, int NrOdfDirections > +OdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> +::OdfMaximaExtractionFilter() + : m_NormalizationMethod(MAX_VEC_NORM) + , m_MaxNumPeaks(2) + , m_RelativePeakThreshold(0.4) + , m_AbsolutePeakThreshold(0) + , m_AngularThreshold(0.9) + , m_NumCoeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder) + , m_Toolkit(MRTRIX) + , m_FlipX(false) + , m_FlipY(false) + , m_FlipZ(false) + , m_ApplyDirectionMatrix(false) + , m_ScaleByGfa(false) + , m_Iterations(10) { - if (m_ShCoeffImage.IsNotNull()) - { - cout << "Using preset coefficient image\n"; - return true; - } - - cout << "Starting qball reconstruction\n"; - try { - QballReconstructionFilterType::Pointer filter = QballReconstructionFilterType::New(); - filter->SetBValue(m_Bvalue); - filter->SetGradientImage( m_DiffusionGradients, m_DiffusionImage ); - filter->SetLambda(0.006); - filter->SetNormalizationMethod(QballReconstructionFilterType::QBAR_SOLID_ANGLE); - filter->Update(); - m_ShCoeffImage = filter->GetCoefficientImage(); - if (m_ShCoeffImage.IsNull()) - return false; - return true; - } - catch (...) - { - return false; - } + this->SetNumberOfRequiredInputs(1); } -// solve ax³ + bx² + cx + d = 0 using cardanos method -template< class TOdfPixelType > -std::vector OdfMaximaExtractionFilter< TOdfPixelType > -::SolveCubic(const double& a, const double& b, const double& c, const double& d) +template< class PixelType, int ShOrder, int NrOdfDirections > +double OdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> +::FindCandidatePeaks(OdfType& odf, double thr, std::vector< DirectionType >& container) { - double A, B, p, q, r, D, offset, ee, tmp, root; - std::vector roots; - double inv3 = 1.0/3.0; - - if (a!=0) // solve ax³ + bx² + cx + d = 0 - { - p = b/a; q = c/a; r = d/a; // x³ + px² + qx + r = 0 - A = q-p*p*inv3; - B = (2.0*p*p*p-9.0*p*q+27.0*r)/27.0; - A = A*inv3; - B = B*0.5; - D = B*B+A*A*A; - offset = p*inv3; - - if (D>0.0) // one real root - { - ee = sqrt(D); - tmp = -B+ee; root = cbrt(tmp); - tmp = -B-ee; root += cbrt(tmp); - root -= offset; roots.push_back(root); - } - else if (D<0.0) // three real roots - { - ee = sqrt(-D); - double tmp2 = -B; - double angle = 2.0*inv3*atan(ee/(sqrt(tmp2*tmp2+ee*ee)+tmp2)); - double sqrt3 = sqrt(3.0); - tmp = cos(angle); - tmp2 = sin(angle); - ee = sqrt(-A); - root = 2*ee*tmp-offset; roots.push_back(root); - root = -ee*(tmp+sqrt3*tmp2)-offset; roots.push_back(root); - root = -ee*(tmp-sqrt3*tmp2)-offset; roots.push_back(root); - } - else // one or two real roots - { - tmp=-B; - tmp=cbrt(tmp); - root=2*tmp-offset; roots.push_back(root); - // TODO: check if this is correct (see history) - if (A!=0 || B!=0) - { - root=-tmp-offset; roots.push_back(root); - } - } - } - else if (b!=0) // solve bx² + cx + d = 0 + double gfa = 1.0; + if (m_ScaleByGfa) + gfa = odf.GetGeneralizedFractionalAnisotropy(); + + //Find the peaks using a finite difference method + bool flag = true; + vnl_vector_fixed< bool, NrOdfDirections > used; used.fill(false); + //Find the peaks + for (int i=0; ithr && gfa*val>m_AbsolutePeakThreshold) // limit to one hemisphere ??? { - D = c*c-4*b*d; - if (D>0) + flag = true; + std::vector< int > neighbours = odf.GetNeighbors(i); + for (unsigned int j=0; j -double OdfMaximaExtractionFilter< TOdfPixelType > -::ODF_dtheta(const double& sn, const double& cs, const double& A, const double& B, const double& C, const double& D, const double& E, const double& F, const double& G, const double& H) +template< class PixelType, int ShOrder, int NrOdfDirections > +void OdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> +::BeforeThreadedGenerateData() { - double dtheta=(G-7*E)*sn*sn + (7*F-35*D-H)*sn*cs + (H+C-F-3*A-5*D)*sn + (0.5*E+B+0.5*G)*cs -0.5*G+3.5*E; - return dtheta; + typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); + itk::Vector spacing = ShCoeffImage->GetSpacing(); + double minSpacing = spacing[0]; + if (spacing[1]GetOrigin(); + itk::Matrix direction = ShCoeffImage->GetDirection(); + ImageRegion<3> imageRegion = ShCoeffImage->GetLargestPossibleRegion(); + + if (m_MaskImage.IsNotNull()) + { + origin = m_MaskImage->GetOrigin(); + direction = m_MaskImage->GetDirection(); + imageRegion = m_MaskImage->GetLargestPossibleRegion(); + } + + itk::Vector spacing3 = ShCoeffImage->GetSpacing(); + itk::Point origin3 = ShCoeffImage->GetOrigin(); + itk::Matrix direction3 = ShCoeffImage->GetDirection(); + itk::ImageRegion<3> imageRegion3 = ShCoeffImage->GetLargestPossibleRegion(); + + itk::Vector spacing4; + itk::Point origin4; + itk::Matrix direction4; + itk::ImageRegion<4> imageRegion4; + + spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; + origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin4[3] = 0; + for (int r=0; r<3; r++) + for (int c=0; c<3; c++) + direction4[r][c] = direction3[r][c]; + direction4[3][3] = 1; + imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); + imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); + imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); + imageRegion4.SetSize(3, m_MaxNumPeaks*3); + + m_PeakImage = PeakImageType::New(); + m_PeakImage->SetSpacing( spacing4 ); + m_PeakImage->SetOrigin( origin4 ); + m_PeakImage->SetDirection( direction4 ); + m_PeakImage->SetRegions( imageRegion4 ); + m_PeakImage->Allocate(); + m_PeakImage->FillBuffer(0.0); + + if (m_MaskImage.IsNull()) + { + m_MaskImage = ItkUcharImgType::New(); + m_MaskImage->SetSpacing( spacing ); + m_MaskImage->SetOrigin( origin ); + m_MaskImage->SetDirection( direction ); + m_MaskImage->SetRegions( imageRegion ); + m_MaskImage->Allocate(); + m_MaskImage->FillBuffer(1); + } + m_NumDirectionsImage = ItkUcharImgType::New(); + m_NumDirectionsImage->SetSpacing( spacing ); + m_NumDirectionsImage->SetOrigin( origin ); + m_NumDirectionsImage->SetDirection( direction ); + m_NumDirectionsImage->SetRegions( imageRegion ); + m_NumDirectionsImage->Allocate(); + m_NumDirectionsImage->FillBuffer(0); + + // calculate SH basis + OdfType odf; + vnl_matrix< double > dir_matrix; + dir_matrix.set_size(3, NrOdfDirections); + for (int i=0; iSetNumberOfThreads(1); } -template< class TOdfPixelType > -double OdfMaximaExtractionFilter< TOdfPixelType > -::ODF_dtheta2(const double& sn, const double& cs, const double& A, const double& B, const double& C, const double& D, const double& E, const double& F, const double& G, const double& H) +template< class PixelType, int ShOrder, int NrOdfDirections > +void OdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> +::AfterThreadedGenerateData() { - double dtheta2=4*(G-7*E)*sn*cs + 2*(7*F-35*D-H)*(2*cs*cs-1) + 2*(H+C-F-3*A-5*D)*cs -(E+2*B+G)*sn; - return dtheta2; -} -template< class TOdfPixelType > -double OdfMaximaExtractionFilter< TOdfPixelType > -::ODF_dphi2(const double& sn, const double& cs, const double& A, const double& B, const double& C, const double& D, const double& E, const double& F, const double& G, const double& H) -{ - double dphi2=35*D*((1+cs)*(1+cs)/4)+(3*A-30*D)*(1+cs)/2.0+3*D-A + 0.5*(7*E*(1+cs)/2.0-3*E+B)*sn + (7*F*(1+cs)/2+C-F)*(1-cs)/2.0 + G*sn*(1-cs)/4.0 + H*((1-cs)*(1-cs)/4); - return dphi2; } -template< class TOdfPixelType > -void OdfMaximaExtractionFilter< TOdfPixelType > -::FindCandidatePeaks(const CoefficientPixelType& SHcoeff) +template< class PixelType, int ShOrder, int NrOdfDirections > +void OdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> +::ThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, ThreadIdType threadID ) { - const double thr = 0.03; // threshold on the derivative of the ODF with respect to theta - const double phi_step = 0.005; // step size for 1D exhaustive search on phi - bool highRes; // when close to maxima increase resolution - double mag, Y, Yp, sn, cs; - double phi, dPhi; - double A, B, C, D, E, F, G, H, Bp, Cp, Ep, Fp, Gp, Hp, Bs, Cs, Es, Fs, Gs, Hs; - CoefficientPixelType a, ap; - a = SHcoeff; ap = SHcoeff; - - m_CandidatePeaks.clear(); // clear peaks of last voxel - - for (int adaptiveStepwidth=0; adaptiveStepwidth<=1; adaptiveStepwidth++) - { - phi=0; - while (phi<(2*itk::Math::pi)) // phi exhaustive search 0..pi - { - // calculate 4th order SH representtaion of ODF and according derivative - for (int l=0; l<=4; l=l+2) - { - for (int m=-l; m<=l; m++) - { - int j=l*(l+1)/2+m; - if (m<0) - { - mag = sqrt(((2*l+1)/(2*itk::Math::pi))*factorial(l+m)/factorial(l-m)); - Y = mag*cos(m*phi); - Yp = -m*mag*sin(m*phi); - } - else if (m==0) - { - Y = sqrt((2*l+1)/(4*itk::Math::pi)); - Yp = 0; - } - else - { - mag = pow(-1.0,m)*sqrt(((2*l+1)/(2*itk::Math::pi))*factorial(l-m)/factorial(l+m)); - Y = mag*sin(m*phi); - Yp = m*mag*cos(m*phi); - } - a[j] = SHcoeff[j]*Y; - ap[j] = SHcoeff[j]*Yp; - } - } - - // ODF - A=0.5*a[3]; B=-3*(a[2]+a[4]); C=3*(a[1]+a[5]); D=0.125*a[10]; E=-2.5*(a[9]+a[11]); - F=7.5*(a[8]+a[12]); G=-105*(a[7]+a[13]); H=105*(a[6]+a[14]); - - // phi derivative - Bp=-3*(ap[2]+ap[4]); Cp=3*(ap[1]+ap[5]); Ep=-2.5*(ap[9]+ap[11]); - Fp=7.5*(ap[8]+ap[12]); Gp=-105*(ap[7]+ap[13]); Hp=105*(ap[6]+ap[14]); - - // 2phi derivative - Bs=-B; Cs=-4*C; Es=-E; - Fs=-4*F; Gs=-9*G; Hs=-16*H; + typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); - // solve cubic for tan(theta) - std::vector tanTheta = SolveCubic(Hp+Cp-Fp, Gp+Bp-3*Ep, 6*Fp+Cp, Bp+4*Ep); + ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, outputRegionForThread ); - highRes = false; - dPhi = phi_step; - - //for each real cubic solution for tan(theta) - for (int n=0; n hessian; - hessian(0,0) = ODF_dtheta2(sn, cs, A, B, C, D, E, F, G, H); - hessian(0,1) = ODF_dtheta(sn, cs, 0, Bp, Cp, 0, Ep, Fp, Gp, Hp); - hessian(1,0) = hessian(0,1); - hessian(1,1) = ODF_dphi2(sn, cs, 0, Bs, Cs, 0, Es, Fs, Gs, Hs); - - double det = vnl_det(hessian); // determinant - double tr = vnl_trace(hessian); // trace - - highRes = true; // we are close to a maximum, so turn on high resolution 1D exhaustive search - if (det>=0 && tr<=0) // check if we really have a local maximum - { - vnl_vector_fixed< double, 2 > peak; - peak[0] = theta; - peak[1] = phi; - m_CandidatePeaks.push_back(peak); - } - } - - if (adaptiveStepwidth) // calculate adaptive step width - { - double t2=tanTheta[n]*tanTheta[n]; double t3=t2*tanTheta[n]; double t4=t3*tanTheta[n]; - double const_step=phi_step*(1+t2)/sqrt(t2+t4+pow((((Hs+Cs-Fs)*t3+(Gs+Bs-3*Es)*t2+(6*Fs+Cs)*tanTheta[n]+(Bs+4*Es))/(3*(Hp+Cp-Fp)*t2+2*(Gp+Bp-3*Ep)*tanTheta[n]+(6*Fp+Cp))),2.0)); - if (const_stepGetPixel(idx3)==0) + { + ++cit; + continue; } -} - -template< class TOdfPixelType > -std::vector< vnl_vector_fixed< double, 3 > > OdfMaximaExtractionFilter< TOdfPixelType > -::ClusterPeaks(const CoefficientPixelType& shCoeff) -{ - const double distThres = 0.4; - int npeaks = 0, nMin = 0; - double dMin, dPos, dNeg, d; - Vector3D u; - std::vector< Vector3D > v; - // initialize container for vector clusters - std::vector < std::vector< Vector3D > > clusters; - clusters.resize(m_CandidatePeaks.size()); + CoefficientPixelType c = cit.Get(); + vnl_vector coeffs; coeffs.set_size((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder); + for (int j=0; j::max(); - for (int n=0; nmax) + max = odf[i]; + } + if (max<0.0001) + { + ++cit; + continue; } - // calculate mean vector of each cluster - for (int i=0; i candidates, final_peaks; + double scale = FindCandidatePeaks(odf, max*m_RelativePeakThreshold*0.9, candidates); // find all local maxima - if (npeaks!=0) + max = 0; + for (unsigned int i=0; i shBasis, sphCoords; - - Cart2Sph(v, sphCoords); // convert candidate peaks to spherical angles - shBasis = CalcShBasis(sphCoords, 4); // evaluate spherical harmonics at each peak - vnl_vector odfVals(npeaks); - odfVals.fill(0.0); - double maxVal = itk::NumericTraits::NonpositiveMin(); - int maxPos; - for (int i=0; imaxVal ) - { - maxVal = odfVals(i); - maxPos = i; - } - } - v.clear(); - std::vector< double > restVals; - for (int i=0; i=m_PeakThreshold*maxVal ) - { - u[0] = odfVals(i)*cos(sphCoords(i,1))*sin(sphCoords(i,0)); - u[1] = odfVals(i)*sin(sphCoords(i,1))*sin(sphCoords(i,0)); - u[2] = odfVals(i)*cos(sphCoords(i,0)); - restVals.push_back(odfVals(i)); - v.push_back(u); - } - npeaks = v.size(); - - if (npeaks>m_MaxNumPeaks) // if still too many peaks, keep only the m_MaxNumPeaks with maximum value - { - std::vector< Vector3D > v2; - for (int i=0; i::NonpositiveMin(); //Get the maximum ODF peak value and the corresponding peak index - for (int i=0; imaxVal ) - { - maxVal = restVals[i]; - maxPos = i; - } - - v2.push_back(v[maxPos]); - restVals[maxPos] = 0; //zero that entry in order to find the next maximum - } - return v2; - } + double spherical[3]; + mitk::sh::Cart2Sph(candidates[i][0], candidates[i][1], candidates[i][2], spherical); + vnl_vector x; + x.set_size(2); + x[0] = spherical[1]; // theta + x[1] = spherical[0]; // phi + + VnlCostFunction cost; + if (m_Toolkit==Toolkit::MRTRIX) + cost.SetProblem(coeffs, ShOrder, true, max); + else + cost.SetProblem(coeffs, ShOrder, false, max); + + vnl_lbfgsb minimizer(cost); + minimizer.set_f_tolerance(1e-6); +// minimizer.set_trace(true); + + vnl_vector l; l.set_size(2); l[0] = 0; l[1] = -itk::Math::pi; + vnl_vector u; u.set_size(2); u[0] = itk::Math::pi; u[1] = itk::Math::pi; + vnl_vector bound_selection; bound_selection.set_size(2); bound_selection.fill(2); + minimizer.set_bound_selection(bound_selection); + minimizer.set_lower_bound(l); + minimizer.set_upper_bound(u); + if (m_Iterations>0) + minimizer.set_max_function_evals(m_Iterations); + minimizer.minimize(x); + + float v = -minimizer.get_end_error()*scale; + candidates[i] = mitk::sh::Sph2Cart(x[0], x[1], v); + if (v>max) + max = v; } - return v; -} -// convert cartesian to spherical coordinates -template< class TOdfPixelType > -void OdfMaximaExtractionFilter< TOdfPixelType > -::Cart2Sph(const std::vector< Vector3D >& dir, vnl_matrix& sphCoords) -{ - sphCoords.set_size(dir.size(), 2); + std::sort( candidates.begin(), candidates.end(), CompareVectors ); // sort peaks - for (int i=0; im_AngularThreshold && val -vnl_matrix OdfMaximaExtractionFilter< TOdfPixelType > -::CalcShBasis(vnl_matrix& sphCoords, const int& shOrder) -{ - int R = (shOrder+1)*(shOrder+2)/2; - int M = sphCoords.rows(); - int j, m; double mag, plm; - vnl_matrix shBasis; - shBasis.set_size(M,R); - for (int p=0; p(l,abs(m),cos(sphCoords(p,0))); - mag = sqrt((double)(2*l+1)/(4.0*itk::Math::pi)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; - - if (m<0) - shBasis(p,j) = sqrt(2.0)*mag*cos(fabs((double)m)*sphCoords(p,1)); - else if (m==0) - shBasis(p,j) = mag; - else - shBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(p,1)); - j++; - } + if (flag) + final_peaks.push_back(v1); } - return shBasis; -} - -template< class TOdfPixelType > -void OdfMaximaExtractionFilter< TOdfPixelType > -::GenerateData() -{ - if (!ReconstructQballImage()) - return; - std::cout << "Starting maxima extraction\n"; + itk::Index<4> idx4; idx4[0] = idx3[0]; idx4[1] = idx3[1]; idx4[2] = idx3[2]; - switch (m_NormalizationMethod) + // fill output image + unsigned int num = final_peaks.size(); + if ( num>m_MaxNumPeaks ) + num = m_MaxNumPeaks; + for (unsigned int i=0; iGetDirection()*dir; + + if (m_FlipX) + dir[0] = -dir[0]; + if (m_FlipY) + dir[1] = -dir[1]; + if (m_FlipZ) + dir[2] = -dir[2]; + + for (unsigned int j = 0; j<3; j++) + { + idx4[3] = i*3 + j; + m_PeakImage->SetPixel(idx4, dir[j]); + } } - - typedef ImageRegionConstIterator< CoefficientImageType > InputIteratorType; - - - InputIteratorType git(m_ShCoeffImage, m_ShCoeffImage->GetLargestPossibleRegion() ); - - itk::Vector spacing = m_ShCoeffImage->GetSpacing(); - double minSpacing = spacing[0]; - if (spacing[1]GetOrigin(); - itk::Matrix direction = m_ShCoeffImage->GetDirection(); - ImageRegion<3> imageRegion = m_ShCoeffImage->GetLargestPossibleRegion(); - - // initialize num directions image - m_NumDirectionsImage = ItkUcharImgType::New(); - m_NumDirectionsImage->SetSpacing( spacing ); - m_NumDirectionsImage->SetOrigin( origin ); - m_NumDirectionsImage->SetDirection( direction ); - m_NumDirectionsImage->SetRegions( imageRegion ); - m_NumDirectionsImage->Allocate(); - m_NumDirectionsImage->FillBuffer(0); - - vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); - vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); - - m_DirectionImageContainer = ItkDirectionImageContainer::New(); - for (int i=0; i nullVec; nullVec.Fill(0.0); - ItkDirectionImage::Pointer img = ItkDirectionImage::New(); - img->SetSpacing( spacing ); - img->SetOrigin( origin ); - img->SetDirection( direction ); - img->SetRegions( imageRegion ); - img->Allocate(); - img->FillBuffer(nullVec); - m_DirectionImageContainer->InsertElement(m_DirectionImageContainer->Size(), img); - } - - if (m_MaskImage.IsNull()) - { - m_MaskImage = ItkUcharImgType::New(); - m_MaskImage->SetSpacing( spacing ); - m_MaskImage->SetOrigin( origin ); - m_MaskImage->SetDirection( direction ); - m_MaskImage->SetRegions( imageRegion ); - m_MaskImage->Allocate(); - m_MaskImage->FillBuffer(1); - } - - itk::ImageRegionIterator dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion()); - itk::ImageRegionIterator maskIt(m_MaskImage, m_MaskImage->GetLargestPossibleRegion()); - - int maxProgress = m_MaskImage->GetLargestPossibleRegion().GetSize()[0]*m_MaskImage->GetLargestPossibleRegion().GetSize()[1]*m_MaskImage->GetLargestPossibleRegion().GetSize()[2]; - - boost::progress_display disp(maxProgress); - - git.GoToBegin(); - while( !git.IsAtEnd() ) - { - ++disp; - if (maskIt.Value()<=0) - { - ++git; - ++dirIt; - ++maskIt; - continue; - } - - CoefficientPixelType c = git.Get(); - FindCandidatePeaks(c); - std::vector< Vector3D > directions = ClusterPeaks(c); - - typename CoefficientImageType::IndexType index = git.GetIndex(); - - float max = 0.0; - for (int i=0; imax) - max = directions.at(i).magnitude(); - if (max<0.0001) - max = 1.0; - - for (int i=0; iGetElement(i); - itk::Vector< float, 3 > pixel; - vnl_vector dir = directions.at(i); - - vtkSmartPointer container = vtkSmartPointer::New(); - itk::ContinuousIndex center; - center[0] = index[0]; - center[1] = index[1]; - center[2] = index[2]; - itk::Point worldCenter; - m_ShCoeffImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); - - switch (m_NormalizationMethod) - { - case NO_NORM: - break; - case SINGLE_VEC_NORM: - dir.normalize(); - break; - case MAX_VEC_NORM: - dir /= max; - break; - } - - dir = m_MaskImage->GetDirection()*dir; - pixel.SetElement(0, dir[0]); - pixel.SetElement(1, dir[1]); - pixel.SetElement(2, dir[2]); - img->SetPixel(index, pixel); - - itk::Point worldStart; - worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing; - worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing; - worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing; - vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); - container->GetPointIds()->InsertNextId(id); - itk::Point worldEnd; - worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing; - worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing; - worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing; - id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); - container->GetPointIds()->InsertNextId(id); - m_VtkCellArray->InsertNextCell(container); - } - - dirIt.Set(directions.size()); - - ++git; - ++dirIt; - ++maskIt; - } - - vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); - directionsPolyData->SetPoints(m_VtkPoints); - directionsPolyData->SetLines(m_VtkCellArray); - m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData); - std::cout << "Maxima extraction finished\n"; + m_NumDirectionsImage->SetPixel(idx3, num); + ++cit; + } + MITK_INFO << "Thread " << threadID << " finished extraction"; } + } #endif // __itkOdfMaximaExtractionFilter_cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.h index 454a19292e..e82badce66 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkOdfMaximaExtractionFilter.h @@ -1,145 +1,182 @@ /*========================================================================= Program: Insight Segmentation & Registration Toolkit Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $ Language: C++ Date: $Date: 2006-03-27 17:01:06 $ Version: $Revision: 1.12 $ Copyright (c) Insight Software Consortium. All rights reserved. See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef __itkOdfMaximaExtractionFilter_h_ #define __itkOdfMaximaExtractionFilter_h_ +#include "itkImageToImageFilter.h" +#include "vnl/vnl_vector_fixed.h" +#include "vnl/vnl_matrix.h" +#include "vnl/algo/vnl_svd.h" +#include "itkVectorContainer.h" +#include "itkVectorImage.h" #include - -#include -#include +#include +#include #include -#include - -namespace itk{ -/** \class OdfMaximaExtractionFilter -Class that estimates the maxima of the 4th order SH representation of an ODF using analytic calculations (according to Aganj et al, MICCAI, 2010) - */ +#include +#include -template< class TOdfPixelType > -class OdfMaximaExtractionFilter : public ProcessObject +class VnlCostFunction : public vnl_cost_function { - public: - enum NormalizationMethods { - NO_NORM, ///< don't normalize peaks - SINGLE_VEC_NORM, ///< normalize peaks to length 1 - MAX_VEC_NORM ///< largest peak is normalized to length 1, other peaks relative to it - }; - - typedef OdfMaximaExtractionFilter Self; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - typedef ProcessObject Superclass; - - /** Method for creation through the object factory. */ - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - - /** Runtime information support. */ - itkTypeMacro(OdfMaximaExtractionFilter, ImageToImageFilter) - - typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballReconstructionFilterType; - typedef QballReconstructionFilterType::CoefficientImageType CoefficientImageType; - typedef CoefficientImageType::PixelType CoefficientPixelType; - typedef itk::VectorImage< short, 3 > DiffusionImageType; - typedef vnl_vector_fixed< double, 3 > Vector3D; - typedef VectorContainer< unsigned int, Vector3D > DirectionContainerType; - typedef VectorContainer< unsigned int, DirectionContainerType::Pointer > ContainerType; - typedef itk::Image ItkUcharImgType; - typedef itk::Image< itk::Vector< float, 3 >, 3> ItkDirectionImage; - typedef itk::VectorContainer< unsigned int, ItkDirectionImage::Pointer > ItkDirectionImageContainer; - - // output - itkGetMacro( OutputFiberBundle, mitk::FiberBundle::Pointer) ///< vector field (peak sizes rescaled for visualization purposes) - itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) ///< number of peaks per voxel - itkGetMacro( DirectionImageContainer, ItkDirectionImageContainer::Pointer) ///< container for output peaks - - // input - itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< only voxels inside the binary mask are processed - itkSetMacro( NormalizationMethod, NormalizationMethods) ///< normalization method of ODF peaks - itkSetMacro( DiffusionGradients, DirectionContainerType::Pointer) ///< input for qball reconstruction - itkSetMacro( DiffusionImage, DiffusionImageType::Pointer) ///< input for qball reconstruction - itkSetMacro( Bvalue, float) ///< input for qball reconstruction - itkSetMacro( ShCoeffImage, CoefficientImageType::Pointer) ///< conatins spherical harmonic coefficients - itkSetMacro( MaxNumPeaks, unsigned int) ///< if more peaks are found, only the largest are kept - itkSetMacro( PeakThreshold, double) ///< threshold on peak length relative to the largest peak in the current voxel - - void GenerateData() override; - -protected: - OdfMaximaExtractionFilter(); - ~OdfMaximaExtractionFilter(){} - - /** CSA Qball reconstruction (SH order 4) **/ - bool ReconstructQballImage(); - - /** calculate roots of cubic equation ax³ + bx² + cx + d = 0 using cardanos method **/ - std::vector SolveCubic(const double& a, const double& b, const double& c, const double& d); - - /** derivatives of SH representation of the ODF **/ - double ODF_dtheta2(const double& sn, const double& cs, const double& A, const double& B, const double& C, const double& D, const double& E, const double& F, const double& G, const double& H); - - /** derivatives of SH representation of the ODF **/ - double ODF_dphi2(const double& sn, const double& cs, const double& A, const double& B, const double& C, const double& D, const double& E, const double& F, const double& G, const double& H); - - /** derivatives of SH representation of the ODF **/ - double ODF_dtheta(const double& sn, const double& cs, const double& A, const double& B, const double& C, const double& D, const double& E, const double& F, const double& G, const double& H); - - /** calculate all directions fullfilling the maximum consitions **/ - void FindCandidatePeaks(const CoefficientPixelType& SHcoeff); - - /** cluster the peaks detected by FindCandidatePeaks and retain maximum m_MaxNumPeaks **/ - std::vector< Vector3D > ClusterPeaks(const CoefficientPixelType& shCoeff); - - void Cart2Sph(const std::vector< Vector3D >& dir, vnl_matrix& sphCoords); - vnl_matrix CalcShBasis(vnl_matrix& sphCoords, const int& shOrder); - - // diffusion weighted image (mandatory input) - DirectionContainerType::Pointer m_DiffusionGradients; ///< input for qball reconstruction - DiffusionImageType::Pointer m_DiffusionImage; ///< input for qball reconstruction - float m_Bvalue; ///< input for qball reconstruction - - // binary mask image (optional input) - ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed + bool mrtrix; + int ShOrder; + vnl_vector coeffs; + double max; + void SetProblem(vnl_vector& coeffs, int ShOrder, bool mrtrix, double max) + { + this->coeffs = coeffs; + this->ShOrder = ShOrder; + this->mrtrix = mrtrix; + this->max = max; + } + + VnlCostFunction(const int NumVars=2) : vnl_cost_function(NumVars) + { + } + + // cost function + double f(vnl_vector const &x) + { + return -mitk::sh::GetValue(coeffs, ShOrder, x[0], x[1], mrtrix); + } + + // gradient of cost function + void gradf(vnl_vector const &x, vnl_vector &dx) + { + fdgradf(x, dx); + } +}; - // input parameters - NormalizationMethods m_NormalizationMethod; ///< normalization for peaks - double m_PeakThreshold; ///< threshold on peak length relative to the largest peak in the current voxel - unsigned int m_MaxNumPeaks; ///< if more peaks are found, only the largest are kept - // intermediate results - CoefficientImageType::Pointer m_ShCoeffImage; ///< conatins spherical harmonic coefficients - std::vector< vnl_vector_fixed< double, 2 > > m_CandidatePeaks; ///< container for candidate peaks (all extrema, also minima) +namespace itk{ - // output data - mitk::FiberBundle::Pointer m_OutputFiberBundle; ///< vector field (peak sizes rescaled for visualization purposes) - ItkUcharImgType::Pointer m_NumDirectionsImage; ///< number of peaks per voxel - ItkDirectionImageContainer::Pointer m_DirectionImageContainer; ///< output peaks +/** +* \brief Extract ODF peaks by searching all local maxima on a roughly sampled sphere and a successive gradient descent optimization +*/ -private: +template< class PixelType, int ShOrder, int NrOdfDirections > +class OdfMaximaExtractionFilter : public ImageToImageFilter< Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, +Image< unsigned char, 3 > > +{ + public: + + enum Toolkit { ///< SH coefficient convention (depends on toolkit) + FSL, + MRTRIX + }; + + enum NormalizationMethods { + NO_NORM, ///< no length normalization of the output peaks + SINGLE_VEC_NORM, ///< normalize the single peaks to length 1 + MAX_VEC_NORM ///< normalize all peaks according to their length in comparison to the largest peak in the respective voxel (0-1) + }; + + typedef OdfMaximaExtractionFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< Image< Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, + Image< unsigned char, 3 > > Superclass; + + /** Method for creation through the object factory. */ + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(OdfMaximaExtractionFilter, ImageToImageFilter) + + typedef typename Superclass::InputImageType CoefficientImageType; + typedef typename CoefficientImageType::RegionType CoefficientImageRegionType; + typedef typename CoefficientImageType::PixelType CoefficientPixelType; + typedef typename Superclass::OutputImageType OutputImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + typedef typename Superclass::InputImageRegionType InputImageRegionType; + typedef mitk::PeakImage::ItkPeakImageType PeakImageType; + + typedef OrientationDistributionFunction OdfType; + typedef itk::Image ItkUcharImgType; + + typedef vnl_vector_fixed< double, 3 > DirectionType; + + // input + itkSetMacro( MaxNumPeaks, unsigned int) ///< maximum number of peaks per voxel. if more peaks are detected, only the largest are kept. + itkSetMacro( RelativePeakThreshold, double) ///< threshold on the peak length relative to the largest peak inside the current voxel + itkSetMacro( AbsolutePeakThreshold, double) ///< hard threshold on the peak length of all local maxima + itkSetMacro( AngularThreshold, double) ///< directions closer together than the specified threshold are discarded + itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< only voxels inside the binary mask are processed + itkSetMacro( NormalizationMethod, NormalizationMethods) ///< normalization method of ODF peaks + itkSetMacro( FlipX, bool) ///< flip peaks in x direction + itkSetMacro( FlipY, bool) ///< flip peaks in y direction + itkSetMacro( FlipZ, bool) ///< flip peaks in z direction + itkSetMacro( ApplyDirectionMatrix, bool) + itkSetMacro( ScaleByGfa, bool) + itkSetMacro( Iterations, int) + + // output + itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer ) + itkGetMacro( PeakImage, PeakImageType::Pointer ) + + itkSetMacro( Toolkit, Toolkit) ///< define SH coefficient convention (depends on toolkit) + itkGetMacro( Toolkit, Toolkit) ///< SH coefficient convention (depends on toolkit) + + protected: + OdfMaximaExtractionFilter(); + ~OdfMaximaExtractionFilter(){} + + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadID ); + void AfterThreadedGenerateData(); + + /** Extract all local maxima from the densely sampled ODF surface. Thresholding possible. **/ + double FindCandidatePeaks(OdfType& odf, double odfMax, std::vector< DirectionType >& inDirs); + + /** Cluster input directions within a certain angular threshold **/ + std::vector< DirectionType > MeanShiftClustering(std::vector< DirectionType >& inDirs); + + private: + + NormalizationMethods m_NormalizationMethod; ///< normalization method of ODF peaks + unsigned int m_MaxNumPeaks; ///< maximum number of peaks per voxel. if more peaks are detected, only the largest are kept. + double m_RelativePeakThreshold; ///< threshold on the peak length relative to the largest peak inside the current voxel + double m_AbsolutePeakThreshold;///< hard threshold on the peak length of all local maxima + vnl_matrix< float > m_ShBasis; ///< container for evaluated SH base functions + double m_AngularThreshold; + const int m_NumCoeffs; ///< number of spherical harmonics coefficients + + PeakImageType::Pointer m_PeakImage; + ItkUcharImgType::Pointer m_NumDirectionsImage; ///< number of peaks per voxel + ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed + + Toolkit m_Toolkit; + bool m_FlipX; + bool m_FlipY; + bool m_FlipZ; + bool m_ApplyDirectionMatrix; + bool m_ScaleByGfa; + int m_Iterations; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkOdfMaximaExtractionFilter.cpp" #endif #endif //__itkOdfMaximaExtractionFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/mitkShImage.h b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/mitkShImage.h index 7d46192a3f..a8a73aaf54 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/mitkShImage.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/mitkShImage.h @@ -1,71 +1,71 @@ /*=================================================================== 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 __mitkShImage__h #define __mitkShImage__h #include "mitkImage.h" #include "itkVectorImage.h" #include "mitkImageVtkAccessor.h" #include namespace mitk { /** * \brief this class encapsulates spherical harmonics coefficients */ class MITKDIFFUSIONCORE_EXPORT ShImage : public Image { public: - typedef itk::Image ShOnDiskType; // we store the sh images in MRtrix 4D float format and convert the on load to 3D multi-component images + typedef itk::Image ShOnDiskType; // we store the sh images in MRtrix 4D float format and convert on load to 3D multi-component images mitkClassMacro( ShImage, Image ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) virtual const vtkImageData* GetNonRgbVtkImageData(int t = 0, int n = 0) const; virtual vtkImageData* GetNonRgbVtkImageData(int t = 0, int n = 0); const vtkImageData* GetVtkImageData(int t = 0, int n = 0) const override; vtkImageData* GetVtkImageData(int t = 0, int n = 0) override; virtual void ConstructRgbImage() const; int ShOrder(); int NumCoefficients(); protected: ShImage(); ~ShImage() override; template void Construct() const; void PrintSelf(std::ostream &os, itk::Indent indent) const override; mutable mitk::Image::Pointer m_RgbImage; int m_ShOrder; int m_NumCoefficients; }; } // namespace mitk #endif /* __mitkShImage__h */ diff --git a/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h b/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h index 22261eac41..2d6a11daf5 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h @@ -1,135 +1,138 @@ /*=================================================================== 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 __mitkDiffusionFunctionCollection_h_ #define __mitkDiffusionFunctionCollection_h_ #include #include #include #include #include #include #include #include namespace mitk{ class MITKDIFFUSIONCORE_EXPORT imv { public: static std::vector< std::pair< itk::Index<3>, double > > IntersectImage(const itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef); static itk::Point GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class TPixelType, class TOutPixelType=TPixelType > static TOutPixelType GetImageValue(const itk::Point& itkP, bool interpolate, typename itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::Pointer interpolator) { if (interpolator==nullptr) return 0.0; itk::ContinuousIndex< float, 3> cIdx; interpolator->ConvertPointToContinuousIndex(itkP, cIdx); if (interpolator->IsInsideBuffer(cIdx)) { if (interpolate) return interpolator->EvaluateAtContinuousIndex(cIdx); else { itk::Index<3> idx; interpolator->ConvertContinuousIndexToNearestIndex(cIdx, idx); return interpolator->EvaluateAtIndex(idx); } } else return 0.0; } template< class TPixelType=unsigned char > static bool IsInsideMask(const itk::Point& itkP, bool interpolate, typename itk::LinearInterpolateImageFunction< itk::Image< TPixelType, 3 >, float >::Pointer interpolator, float threshold=0.5) { if (interpolator==nullptr) return false; itk::ContinuousIndex< float, 3> cIdx; interpolator->ConvertPointToContinuousIndex(itkP, cIdx); if (interpolator->IsInsideBuffer(cIdx)) { double value = 0.0; if (interpolate) value = interpolator->EvaluateAtContinuousIndex(cIdx); else { itk::Index<3> idx; interpolator->ConvertContinuousIndexToNearestIndex(cIdx, idx); value = interpolator->EvaluateAtIndex(idx); } if (value>=threshold) return true; } return false; } }; class MITKDIFFUSIONCORE_EXPORT sh { public: static double factorial(int number); static void Cart2Sph(double x, double y, double z, double* spherical); + static vnl_vector_fixed Sph2Cart(const double& theta, const double& phi, const double& rad); static double legendre0(int l); static double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart); static double Yj(int m, int k, float theta, float phi, bool mrtrix=true); static vnl_matrix CalcShBasisForDirections(int sh_order, vnl_matrix U, bool mrtrix=true); + static float GetValue(const vnl_vector& coefficients, const int& sh_order, const vnl_vector_fixed& dir, const bool mrtrix); + static float GetValue(const vnl_vector &coefficients, const int &sh_order, const double theta, const double phi, const bool mrtrix); }; class MITKDIFFUSIONCORE_EXPORT gradients { private: typedef std::vector IndiciesVector; typedef mitk::BValueMapProperty::BValueMap BValueMap; typedef DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; typedef DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; public: static GradientDirectionContainerType::Pointer ReadBvalsBvecs(std::string bvals_file, std::string bvecs_file, double& reference_bval); static void WriteBvalsBvecs(std::string bvals_file, std::string bvecs_file, GradientDirectionContainerType::Pointer gradients, double reference_bval); static std::vector GetAllUniqueDirections(const BValueMap &bValueMap, GradientDirectionContainerType *refGradientsContainer ); static bool CheckForDifferingShellDirections(const BValueMap &bValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer); static vnl_matrix ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder); static vnl_matrix ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer); static GradientDirectionContainerType::Pointer CreateNormalizedUniqueGradientDirectionContainer(const BValueMap &bValueMap, const GradientDirectionContainerType * origninalGradentcontainer); }; } #endif //__mitkDiffusionFunctionCollection_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp index db526ca271..2279e3df08 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp @@ -1,508 +1,550 @@ /*=================================================================== 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 "mitkDiffusionFunctionCollection.h" #include "mitkNumericTypes.h" #include #include #include #include #include "itkVectorContainer.h" #include "vnl/vnl_vector.h" #include #include #include #include // Intersect a finite line (with end points p0 and p1) with all of the // cells of a vtkImageData std::vector< std::pair< itk::Index<3>, double > > mitk::imv::IntersectImage(const itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef) { std::vector< std::pair< itk::Index<3>, double > > out; if (si == ei) { double d[3]; for (int i=0; i<3; ++i) d[i] = (sf[i]-ef[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); out.push_back( std::pair< itk::Index<3>, double >(si, len) ); return out; } double bounds[6]; double entrancePoint[3]; double exitPoint[3]; double startPoint[3]; double endPoint[3]; double t0, t1; for (int i=0; i<3; ++i) { startPoint[i] = sf[i]; endPoint[i] = ef[i]; if (si[i]>ei[i]) { int t = si[i]; si[i] = ei[i]; ei[i] = t; } } for (int x = si[0]; x<=ei[0]; ++x) for (int y = si[1]; y<=ei[1]; ++y) for (int z = si[2]; z<=ei[2]; ++z) { bounds[0] = (double)x - 0.5; bounds[1] = (double)x + 0.5; bounds[2] = (double)y - 0.5; bounds[3] = (double)y + 0.5; bounds[4] = (double)z - 0.5; bounds[5] = (double)z + 0.5; int entryPlane; int exitPlane; int hit = vtkBox::IntersectWithLine(bounds, startPoint, endPoint, t0, t1, entrancePoint, exitPoint, entryPlane, exitPlane); if (hit) { if (entryPlane>=0 && exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i] - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (entryPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (ef[i] - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i]-sf[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } } } return out; } //------------------------- SH-function ------------------------------------ double mitk::sh::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } void mitk::sh::Cart2Sph(double x, double y, double z, double *spherical) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); if( rad < mitk::eps ) { th = itk::Math::pi/2; phi = itk::Math::pi/2; } else { th = acos(z/rad); phi = atan2(y, x); } spherical[0] = phi; spherical[1] = th; spherical[2] = rad; } +vnl_vector_fixed mitk::sh::Sph2Cart(const double& theta, const double& phi, const double& rad) +{ + vnl_vector_fixed dir; + dir[0] = rad * sin(theta) * cos(phi); + dir[1] = rad * sin(theta) * sin(phi); + dir[2] = rad * cos(theta); + return dir; +} + double mitk::sh::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i(l,abs(m),-cos(theta)); double mag = sqrt((double)(2*l+1)/(4.0*itk::Math::pi)*::boost::math::factorial(l-abs(m))/::boost::math::factorial(l+abs(m)))*plm; if (m>0) return mag*cos(m*phi); else if (m==0) return mag; else return mag*sin(-m*phi); } return 0; } vnl_matrix mitk::sh::CalcShBasisForDirections(int sh_order, vnl_matrix U, bool mrtrix) { vnl_matrix sh_basis = vnl_matrix(U.cols(), (sh_order*sh_order + sh_order + 2)/2 + sh_order ); for(unsigned int i=0; i &coefficients, const int &sh_order, const double theta, const double phi, const bool mrtrix) +{ + float val = 0; + for(int k=0; k<=sh_order; k+=2) + { + for(int m=-k; m<=k; m++) + { + int j = (k*k + k + 2)/2 + m - 1; + val += coefficients[j] * mitk::sh::Yj(m, k, theta, phi, mrtrix); + } + } + + return val; +} + +float mitk::sh::GetValue(const vnl_vector &coefficients, const int &sh_order, const vnl_vector_fixed &dir, const bool mrtrix) +{ + double spherical[3]; + mitk::sh::Cart2Sph(dir[0], dir[1], dir[2], spherical); + + float val = 0; + for(int k=0; k<=sh_order; k+=2) + { + for(int m=-k; m<=k; m++) + { + int j = (k*k + k + 2)/2 + m - 1; + val += coefficients[j] * mitk::sh::Yj(m, k, spherical[1], spherical[0], mrtrix); + } + } + + return val; +} + //------------------------- gradients-function ------------------------------------ mitk::gradients::GradientDirectionContainerType::Pointer mitk::gradients::ReadBvalsBvecs(std::string bvals_file, std::string bvecs_file, double& reference_bval) { mitk::gradients::GradientDirectionContainerType::Pointer directioncontainer = mitk::gradients::GradientDirectionContainerType::New(); std::vector bvec_entries; if (!itksys::SystemTools::FileExists(bvecs_file)) mitkThrow() << "bvecs file not existing: " << bvecs_file; else { std::string line; std::ifstream myfile (bvecs_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; boost::split(strs,line,boost::is_any_of("\t \n")); for (auto token : strs) { if (!token.empty()) { try { bvec_entries.push_back(boost::lexical_cast(token)); } catch(...) { mitkThrow() << "Encountered invalid bvecs file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvecs file could not be opened: " << bvals_file; } } reference_bval = -1; std::vector bval_entries; if (!itksys::SystemTools::FileExists(bvals_file)) mitkThrow() << "bvals file not existing: " << bvals_file; else { std::string line; std::ifstream myfile (bvals_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; boost::split(strs,line,boost::is_any_of("\t \n")); for (auto token : strs) { if (!token.empty()) { try { bval_entries.push_back(boost::lexical_cast(token)); if (bval_entries.back()>reference_bval) reference_bval = bval_entries.back(); } catch(...) { mitkThrow() << "Encountered invalid bvals file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvals file could not be opened: " << bvals_file; } } for(unsigned int i=0; i 0) { vec.normalize(); vec[0] = sqrt(factor)*vec[0]; vec[1] = sqrt(factor)*vec[1]; vec[2] = sqrt(factor)*vec[2]; } directioncontainer->InsertElement(i,vec); } return directioncontainer; } void mitk::gradients::WriteBvalsBvecs(std::string bvals_file, std::string bvecs_file, GradientDirectionContainerType::Pointer gradients, double reference_bval) { std::ofstream myfile; myfile.open (bvals_file.c_str()); for(unsigned int i=0; iSize(); i++) { double twonorm = gradients->ElementAt(i).two_norm(); myfile << std::round(reference_bval*twonorm*twonorm) << " "; } myfile.close(); std::ofstream myfile2; myfile2.open (bvecs_file.c_str()); for(int j=0; j<3; j++) { for(unsigned int i=0; iSize(); i++) { GradientDirectionType direction = gradients->ElementAt(i); direction.normalize(); myfile2 << direction.get(j) << " "; } myfile2 << std::endl; } } std::vector mitk::gradients::GetAllUniqueDirections(const BValueMap & refBValueMap, GradientDirectionContainerType *refGradientsContainer ) { IndiciesVector directioncontainer; auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); auto containerIt = directioncontainer.begin(); bool directionExist = false; while(containerIt != directioncontainer.end()) { if (fabs(dot_product(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { directioncontainer.push_back(wntIndex); } } } return directioncontainer; } bool mitk::gradients::CheckForDifferingShellDirections(const BValueMap & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer) { auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ auto mapIterator_2 = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator_2++; //skip bzero Values for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){ if(mapIterator_2 == mapIterator) continue; IndiciesVector currentShell = mapIterator->second; IndiciesVector testShell = mapIterator_2->second; for (unsigned int i = 0; i< currentShell.size(); i++) if (fabs(dot_product(refGradientsContainer->ElementAt(currentShell[i]), refGradientsContainer->ElementAt(testShell[i]))) <= 0.9998) { return true; } } } return false; } vnl_matrix mitk::gradients::ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer) { vnl_matrix Q(3, refShell.size()); Q.fill(0.0); for(unsigned int i = 0; i < refShell.size(); i++) { GradientDirectionType dir = refGradientsContainer->ElementAt(refShell[i]); double x = dir.normalize().get(0); double y = dir.normalize().get(1); double z = dir.normalize().get(2); double cart[3]; mitk::sh::Cart2Sph(x,y,z,cart); Q(0,i) = cart[0]; Q(1,i) = cart[1]; Q(2,i) = cart[2]; } return Q; } vnl_matrix mitk::gradients::ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder) { vnl_matrix SHBasisOutput(QBallReference.cols(), (LOrder+1)*(LOrder+2)*0.5); SHBasisOutput.fill(0.0); for(int i=0; i< (int)SHBasisOutput.rows(); i++) for(int k = 0; k <= (int)LOrder; k += 2) for(int m =- k; m <= k; m++) { int j = ( k * k + k + 2 ) / 2.0 + m - 1; double phi = QBallReference(0,i); double th = QBallReference(1,i); double val = mitk::sh::Yj(m,k,th,phi); SHBasisOutput(i,j) = val; } return SHBasisOutput; } mitk::gradients::GradientDirectionContainerType::Pointer mitk::gradients::CreateNormalizedUniqueGradientDirectionContainer(const mitk::gradients::BValueMap & bValueMap, const GradientDirectionContainerType *origninalGradentcontainer) { mitk::gradients::GradientDirectionContainerType::Pointer directioncontainer = mitk::gradients::GradientDirectionContainerType::New(); auto mapIterator = bValueMap.begin(); if(bValueMap.find(0) != bValueMap.end() && bValueMap.size() > 1){ mapIterator++; //skip bzero Values vnl_vector_fixed vec; vec.fill(0.0); directioncontainer->push_back(vec); } for( ; mapIterator != bValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); mitk::gradients::GradientDirectionContainerType::Iterator containerIt = directioncontainer->Begin(); bool directionExist = false; while(containerIt != directioncontainer->End()) { if (fabs(dot_product(containerIt.Value(), origninalGradentcontainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { GradientDirectionType dir(origninalGradentcontainer->ElementAt(wntIndex)); directioncontainer->push_back(dir.normalize()); } } } return directioncontainer; } diff --git a/Modules/DiffusionImaging/DiffusionIO/mitkPeakImageMapper2D.cpp b/Modules/DiffusionImaging/DiffusionIO/mitkPeakImageMapper2D.cpp index 1337f34d8f..ca85bf8375 100644 --- a/Modules/DiffusionImaging/DiffusionIO/mitkPeakImageMapper2D.cpp +++ b/Modules/DiffusionImaging/DiffusionIO/mitkPeakImageMapper2D.cpp @@ -1,179 +1,180 @@ /*=================================================================== 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 "mitkPeakImageMapper2D.h" #include "mitkBaseRenderer.h" #include "mitkDataNode.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include mitk::PeakImageMapper2D::PeakImageMapper2D() { m_lut = vtkSmartPointer::New(); m_lut->Build(); } mitk::PeakImageMapper2D::~PeakImageMapper2D() { } mitk::PeakImage* mitk::PeakImageMapper2D::GetInput() { return dynamic_cast< mitk::PeakImage * > ( GetDataNode()->GetData() ); } void mitk::PeakImageMapper2D::UpdateVtkTransform(mitk::BaseRenderer *) { // don't apply transform since the peak polydata is already in world coordinates. return; } void mitk::PeakImageMapper2D::Update(mitk::BaseRenderer * renderer) { mitk::DataNode* node = this->GetDataNode(); if (node == nullptr) return; bool visible = true; node->GetVisibility(visible, renderer, "visible"); if ( !visible ) return; this->GenerateDataForRenderer( renderer ); } // vtkActors and Mappers are feeded here void mitk::PeakImageMapper2D::GenerateDataForRenderer(mitk::BaseRenderer *renderer) { mitk::PeakImage* peakImage = this->GetInput(); //the handler of local storage gets feeded in this method with requested data for related renderwindow LocalStorage *localStorage = m_LocalStorageHandler.GetLocalStorage(renderer); vtkSmartPointer polyData = peakImage->GetPolyData(); if (polyData == nullptr) return; localStorage->m_Mapper->ScalarVisibilityOn(); localStorage->m_Mapper->SetScalarModeToUsePointFieldData(); localStorage->m_Mapper->SetLookupTable(m_lut); //apply the properties after the slice was set localStorage->m_Mapper->SelectColorArray("FIBER_COLORS"); localStorage->m_Mapper->SetInputData(polyData); mitk::SliceNavigationController::Pointer sliceContr = renderer->GetSliceNavigationController(); mitk::PlaneGeometry::ConstPointer planeGeo = sliceContr->GetCurrentPlaneGeometry(); - mitk::Point3D plane_origin = planeGeo->GetCenter(); + mitk::Point3D plane_origin = planeGeo->GetOrigin(); mitk::DataNode* node = this->GetDataNode(); mitk::Image* image = dynamic_cast(node->GetData()); mitk::Vector3D spacing = image->GetGeometry()->GetSpacing(); localStorage->m_Mapper->RemoveAllClippingPlanes(); float clipping_plane_thickness = 1; if(spacing[0]GetNormal(); + plane_normal.Normalize(); double vnormal[3]; double vp1[3]; double vp2[3]; - vp1[0] = plane_origin[0] + plane_normal[0] * clipping_plane_thickness; - vp1[1] = plane_origin[1] + plane_normal[1] * clipping_plane_thickness; - vp1[2] = plane_origin[2] + plane_normal[2] * clipping_plane_thickness; + vp1[0] = plane_origin[0] - plane_normal[0] * clipping_plane_thickness; + vp1[1] = plane_origin[1] - plane_normal[1] * clipping_plane_thickness; + vp1[2] = plane_origin[2] - plane_normal[2] * clipping_plane_thickness; - vp2[0] = plane_origin[0] - plane_normal[0] * clipping_plane_thickness; - vp2[1] = plane_origin[1] - plane_normal[1] * clipping_plane_thickness; - vp2[2] = plane_origin[2] - plane_normal[2] * clipping_plane_thickness; + vp2[0] = plane_origin[0] + plane_normal[0] * clipping_plane_thickness; + vp2[1] = plane_origin[1] + plane_normal[1] * clipping_plane_thickness; + vp2[2] = plane_origin[2] + plane_normal[2] * clipping_plane_thickness; { vnormal[0] = vp2[0] - vp1[0]; vnormal[1] = vp2[1] - vp1[1]; vnormal[2] = vp2[2] - vp1[2]; vtkSmartPointer plane = vtkSmartPointer::New(); plane->SetOrigin(vp1); plane->SetNormal(vnormal); localStorage->m_Mapper->AddClippingPlane(plane); } { vnormal[0] = vp1[0] - vp2[0]; vnormal[1] = vp1[1] - vp2[1]; vnormal[2] = vp1[2] - vp2[2]; vtkSmartPointer plane = vtkSmartPointer::New(); plane->SetOrigin(vp2); plane->SetNormal(vnormal); localStorage->m_Mapper->AddClippingPlane(plane); } localStorage->m_PointActor->SetMapper(localStorage->m_Mapper); float linewidth = 1.0; this->GetDataNode()->GetFloatProperty("shape.linewidth",linewidth); localStorage->m_PointActor->GetProperty()->SetLineWidth(linewidth); // We have been modified => save this for next Update() localStorage->m_LastUpdateTime.Modified(); } vtkProp* mitk::PeakImageMapper2D::GetVtkProp(mitk::BaseRenderer *renderer) { this->Update(renderer); return m_LocalStorageHandler.GetLocalStorage(renderer)->m_PointActor; } void mitk::PeakImageMapper2D::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer, bool overwrite) { Superclass::SetDefaultProperties(node, renderer, overwrite); //add other parameters to propertylist node->AddProperty( "color", mitk::ColorProperty::New(1.0,1.0,1.0), renderer, overwrite); node->AddProperty( "shape.linewidth", mitk::FloatProperty::New(1.0), renderer, overwrite); } mitk::PeakImageMapper2D::LocalStorage::LocalStorage() { m_PointActor = vtkSmartPointer::New(); m_Mapper = vtkSmartPointer::New(); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox index aad2127b12..5cdb039c00 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox @@ -1,44 +1,9 @@ /** \page org_mitk_views_odfmaximaextraction ODF Peak Extraction -This view provides the user interface to extract the principal diffusion direction of tensors and the peaks of spherical harmonic ODFs. It also allows to import SH-coefficient images from MRtrix and FSL to an MITK-readable format. +This view provides the user interface to extract the principal diffusion direction of tensors and the peaks of spherical harmonic ODFs. -The output peaks and the imported ODF images can for example be used for streamline tractography. +The output peaks can for example be used for streamline tractography. -Available sections: - - \ref OdfMaxUserManualInputData - - \ref OdfMaxUserManualOutputData - - \ref OdfMaxUserManualMethods - - \ref OdfMaxUserManualParameters - -\section OdfMaxUserManualInputData Input Data - -Mandatory Input: - -\li Tensor image or image containing spherical harmonic coefficients. The SH coefficient images can be obtain from the Q-Ball reconstruction view. - -Optional Input: - -\li Binary mask to define the extraction area. - -\section OdfMaxUserManualOutputData Output Data - -\li Vector field: 3D representation of the resulting peaks. Only for visualization purposes (the peaks are scaled additionally to the specified normalization to improve the visualization)! -\li Peak Image: 4D image containing the peak information. -\li Num. of Peaks per Voxel: Image containing the number of extracted peaks per voxel as image value. - -\section OdfMaxUserManualMethods Peak Extraction Methods - -\li If a tensor image is used as input, the output is simply the largest eigenvector of each voxel. -\li If a SH coefficient image is used as input, all local maxima of the densely sampled sphere are extracted and resulting directions that point in a similar direction are clustered to obtain the principal directions. - -\section OdfMaxUserManualParameters Input Parameters - -\li Vector normalization method (no normalization, maximum normalization of the vecors of one voxel and independent normalization of each vecor). -\li Maximum number of peaks to extract. If more peaks are found only the largest are kept. -\li Relative threshold to discard small peaks. Value relative to the largest peak of the respective voxel. -\li Absolute threshold on the peak size. To evaluate this threshold the peaks are additionally weighted by their GFA (low GFA voxels are more likely to be discarded). -\li Clustering Angle: close peaks are clustered to obtain the principal direction. -\li Angular Threshold: peaks that were too far away to be clustered can be discarded by this threshold. */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp index 5e9242917c..6990aec704 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp @@ -1,380 +1,380 @@ /*=================================================================== 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. ===================================================================*/ //misc #include // Blueberry #include #include // Qmitk #include "QmitkOdfMaximaExtractionView.h" #include #include #include #include #include #include #include -#include +#include #include #include #include #include #include #include #include #include #include #include #include // Qt #include const std::string QmitkOdfMaximaExtractionView::VIEW_ID = "org.mitk.views.odfmaximaextractionview"; using namespace mitk; QmitkOdfMaximaExtractionView::QmitkOdfMaximaExtractionView() : m_Controls(nullptr) { } // Destructor QmitkOdfMaximaExtractionView::~QmitkOdfMaximaExtractionView() { } void QmitkOdfMaximaExtractionView::CreateQtPartControl(QWidget *parent) { // build up qt view, unless already done if (!m_Controls) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkOdfMaximaExtractionViewControls; m_Controls->setupUi(parent); connect((QObject*)m_Controls->m_StartPeakExtractionButton, SIGNAL(clicked()), (QObject*) this, SLOT(StartPeakExtraction())); connect((QObject*)m_Controls->m_ImportShCoeffs, SIGNAL(clicked()), (QObject*) this, SLOT(ConvertShCoeffs())); m_Controls->m_MaskBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateNot::Pointer isDwi = mitk::NodePredicateNot::New(mitk::NodePredicateIsDWI::New()); mitk::NodePredicateNot::Pointer isOdf = mitk::NodePredicateNot::New(mitk::NodePredicateDataType::New("OdfImage")); mitk::NodePredicateAnd::Pointer unwanted = mitk::NodePredicateAnd::New(isOdf, isDwi); mitk::NodePredicateDimension::Pointer dim3 = mitk::NodePredicateDimension::New(3); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); m_Controls->m_MaskBox->SetPredicate(mitk::NodePredicateAnd::New(mitk::NodePredicateAnd::New(unwanted, dim3), isBinaryPredicate)); m_Controls->m_ImageBox->SetPredicate(mitk::NodePredicateAnd::New(mitk::NodePredicateAnd::New(unwanted, isMitkImage), mitk::NodePredicateNot::New(isBinaryPredicate))); m_Controls->m_MaskBox->SetZeroEntryText("--"); connect( (QObject*)(m_Controls->m_ImageBox), SIGNAL(OnSelectionChanged(const mitk::DataNode*)), this, SLOT(OnImageSelectionChanged()) ); m_Controls->m_StartPeakExtractionButton->setVisible(false); m_Controls->m_ImportShCoeffs->setVisible(false); } } void QmitkOdfMaximaExtractionView::SetFocus() { } void QmitkOdfMaximaExtractionView::StartPeakExtraction() { if (dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()) != nullptr) { StartTensorPeakExtraction(dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData())); } else { StartMaximaExtraction(dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData())); } } template void QmitkOdfMaximaExtractionView::TemplatedConvertShCoeffs(mitk::Image* mitkImg) { typedef itk::ShCoefficientImageImporter< float, shOrder > FilterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImg); caster->Update(); typename FilterType::Pointer filter = FilterType::New(); filter->SetInputImage(caster->GetOutput()); filter->GenerateData(); typename FilterType::CoefficientImageType::Pointer itkCi = filter->GetCoefficientImage(); { mitk::Image::Pointer img = dynamic_cast(mitk::ShImage::New().GetPointer()); img->InitializeByItk(itkCi.GetPointer()); img->SetVolume(itkCi->GetBufferPointer()); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_ShImage_Imported"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); } } void QmitkOdfMaximaExtractionView::ConvertShCoeffs() { if (m_Controls->m_ImageBox->GetSelectedNode().IsNull()) return; Image::Pointer mitkImg = dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); if (mitkImg->GetDimension() != 4 && mitkImg->GetLargestPossibleRegion().GetSize()[3]<6) { MITK_INFO << "wrong image type (need 4 dimensions)"; return; } int nrCoeffs = mitkImg->GetLargestPossibleRegion().GetSize()[3]; switch (nrCoeffs) { case 6: TemplatedConvertShCoeffs<2>(mitkImg); break; case 15: TemplatedConvertShCoeffs<4>(mitkImg); break; case 28: TemplatedConvertShCoeffs<6>(mitkImg); break; case 45: TemplatedConvertShCoeffs<8>(mitkImg); break; case 66: TemplatedConvertShCoeffs<10>(mitkImg); break; case 91: TemplatedConvertShCoeffs<12>(mitkImg); break; default : QMessageBox::warning(nullptr, "Error", "Only spherical harmonics orders 2-12 are supported.", QMessageBox::Ok); } } void QmitkOdfMaximaExtractionView::StartTensorPeakExtraction(mitk::TensorImage* img) { typedef itk::DiffusionTensorPrincipalDirectionImageFilter< float > MaximaExtractionFilterType; MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); filter->SetUsePolarCoordinates(false); mitk::BaseGeometry::Pointer geometry; try{ ItkTensorImage::Pointer itkImage = ItkTensorImage::New(); CastToItkImage(img, itkImage); filter->SetInput(itkImage); geometry = img->GetGeometry(); } catch (itk::ExceptionObject &e) { MITK_INFO << "wrong image type: " << e.what(); QMessageBox::warning(nullptr, "Wrong pixel type", "Could not perform Tensor Principal Direction Extraction due to Image has wrong pixel type.", QMessageBox::Ok); return; } if (m_Controls->m_MaskBox->GetSelectedNode().IsNotNull()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); Image::Pointer mitkMaskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); CastToItkImage(mitkMaskImg, itkMaskImage); filter->SetMaskImage(itkMaskImage); } if (m_Controls->m_NormalizationBox->currentIndex() == 0) filter->SetNormalizeVectors(false); filter->SetFaThreshold(m_Controls->m_AbsoluteThresholdBox->value()); filter->Update(); MaximaExtractionFilterType::PeakImageType::Pointer itkImg = filter->GetPeakImage(); mitk::Image::Pointer mitkPeakImage = dynamic_cast(PeakImage::New().GetPointer()); CastToMitkImage(itkImg, mitkPeakImage); DataNode::Pointer node = DataNode::New(); node->SetData(mitkPeakImage); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_PrincipalDirection"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); if (m_Controls->m_OutputNumDirectionsBox->isChecked()) { ItkUcharImgType::Pointer numDirImage = filter->GetOutput(); mitk::Image::Pointer image2 = mitk::Image::New(); image2->InitializeByItk(numDirImage.GetPointer()); image2->SetVolume(numDirImage->GetBufferPointer()); DataNode::Pointer node2 = DataNode::New(); node2->SetData(image2); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_NumDirections"; node2->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node2, m_Controls->m_ImageBox->GetSelectedNode()); } } template void QmitkOdfMaximaExtractionView::StartMaximaExtraction(Image *image) { - typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 10000 > MaximaExtractionFilterType; + typedef itk::OdfMaximaExtractionFilter< float, shOrder, ODF_SAMPLING_SIZE > MaximaExtractionFilterType; typename MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); switch (m_Controls->m_ToolkitBox->currentIndex()) { case 0: filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); break; case 1: filter->SetToolkit(MaximaExtractionFilterType::FSL); break; default: filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); } mitk::BaseGeometry::Pointer geometry; try{ typedef ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); filter->SetInput(caster->GetOutput()); geometry = image->GetGeometry(); } catch (itk::ExceptionObject &e) { MITK_INFO << "wrong image type: " << e.what(); QMessageBox::warning(nullptr, "Wrong pixel type", "Could not perform Finite Differences Extraction due to Image has wrong pixel type.", QMessageBox::Ok); return; } filter->SetAngularThreshold(cos((float)m_Controls->m_AngularThreshold->value()*itk::Math::pi / 180)); - filter->SetClusteringThreshold(cos((float)m_Controls->m_ClusteringAngleBox->value()*itk::Math::pi / 180)); filter->SetMaxNumPeaks(m_Controls->m_MaxNumPeaksBox->value()); - filter->SetPeakThreshold(m_Controls->m_PeakThresholdBox->value()); + filter->SetRelativePeakThreshold(m_Controls->m_PeakThresholdBox->value()); filter->SetAbsolutePeakThreshold(m_Controls->m_AbsoluteThresholdBox->value()); + filter->SetScaleByGfa(m_Controls->m_ScaleByGfaBox->isChecked()); if (m_Controls->m_MaskBox->GetSelectedNode().IsNotNull()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); Image::Pointer mitkMaskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); CastToItkImage(mitkMaskImg, itkMaskImage); filter->SetMaskImage(itkMaskImage); } switch (m_Controls->m_NormalizationBox->currentIndex()) { case 0: filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } filter->Update(); typename MaximaExtractionFilterType::PeakImageType::Pointer itkImg = filter->GetPeakImage(); mitk::Image::Pointer img = dynamic_cast(PeakImage::New().GetPointer()); CastToMitkImage(itkImg, img); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_PEAKS"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); if (m_Controls->m_OutputNumDirectionsBox->isChecked()) { ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); mitk::Image::Pointer image2 = mitk::Image::New(); CastToMitkImage(numDirImage, image2); DataNode::Pointer node2 = DataNode::New(); node2->SetData(image2); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_NUM_DIRECTIONS"; node2->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node2, m_Controls->m_ImageBox->GetSelectedNode()); } } void QmitkOdfMaximaExtractionView::StartMaximaExtraction(Image* img) { mitk::PixelType pixT = img->GetPixelType(); switch (pixT.GetNumberOfComponents()) { case 6: StartMaximaExtraction<2>(img); break; case 15: StartMaximaExtraction<4>(img); break; case 28: StartMaximaExtraction<6>(img); break; case 45: StartMaximaExtraction<8>(img); break; case 66: StartMaximaExtraction<10>(img); break; case 91: StartMaximaExtraction<12>(img); break; default : QMessageBox::warning(nullptr, "Error", "Only spherical harmonics orders 2-12 are supported.", QMessageBox::Ok); } } void QmitkOdfMaximaExtractionView::OnSelectionChanged(berry::IWorkbenchPart::Pointer , const QList& nodes) { (void) nodes; this->OnImageSelectionChanged(); } void QmitkOdfMaximaExtractionView::OnImageSelectionChanged() { m_Controls->m_StartPeakExtractionButton->setVisible(false); m_Controls->m_ImportShCoeffs->setVisible(false); mitk::DataNode::Pointer node = m_Controls->m_ImageBox->GetSelectedNode(); if (node.IsNull()) return; Image::Pointer img = dynamic_cast(node->GetData()); if (img->GetDimension()==4) m_Controls->m_ImportShCoeffs->setVisible(true); else m_Controls->m_StartPeakExtractionButton->setVisible(true); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionViewControls.ui index 80d9d455b5..f2164ff799 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionViewControls.ui @@ -1,476 +1,470 @@ - + QmitkOdfMaximaExtractionViewControls 0 0 397 848 Form QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 25 QFrame::NoFrame QFrame::Raised 0 0 0 0 true Generate ODF image and MITK compatible SH coefficient from other toolkits. Start SH Coefficient Import true Extract ODF peaks using finite differences on the densely sampled ODF surface. Start Peak Extraction Please Select Input Data 6 6 6 6 Select a tensor image or a SH coefficient image (generate using Q-Ball reconstruction view). ShCoeff/DTI: Mask Image: - + - + Additional Output QFormLayout::AllNonFixedFieldsGrow 6 6 6 6 Output unsigned char image containing the number of directions per voxel. #Peaks per Voxel false Parameters 6 6 6 6 QFrame::NoFrame QFrame::Raised 0 0 0 0 6 Vector Normalization: <html><head/><body><p>The vector fields are always coorected for image spacing and using the lagest eigenvalue in case of the tensor peak extraction. This is done for visualizytion purposes. The output direction images are not affected.</p></body></html> 1 No Normalization MAX Normalize Single Vec Normalization true QFrame::NoFrame QFrame::Raised 0 0 0 0 6 - - - - true - - - Max. Peaks: - - - - - - - Relative Threshold: - - - true Peak threshold relative to the largest peak per voxel. 3 0.000000000000000 1.000000000000000 0.100000000000000 - 0.500000000000000 + 0.400000000000000 + + + + + + + true + + + Max. Peaks: true Absolute peak threshold (only used for the finite differences method). The value is additionally scaled by 1/GFA. 3 0.000000000000000 1.000000000000000 0.010000000000000 0.030000000000000 true Maximum number of peaks to extract. 1 1000 3 - + - Clustering Angle: + Angular Threshold: - + - Cluster close directions. Define "close" here. + Discard smaller peaks in the defined angle around the maximum peaks that were too far away to be clustered. + + + 0 90 - 30 + 15 Absolute Threshold: + + + + Relative Threshold: + + + - + - Angular Threshold: + Scale by GFA: - - - Discard smaller peaks in the defined angle around the maximum peaks that were too far away to be clustered. - - - 0 - - - 90 - - - 0 + + + Qt::Vertical 20 259 Spherical Harmonic Convention 6 6 6 6 Define SH coefficient convention (depends on toolkit) 0 MITK/MRtrix FSL QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 - - - \ No newline at end of file + + +