diff --git a/Modules/DICOMQI/autoload/IO/mitkDICOMSegmentationIO.cpp b/Modules/DICOMQI/autoload/IO/mitkDICOMSegmentationIO.cpp index b269692655..c0b5d43ab3 100644 --- a/Modules/DICOMQI/autoload/IO/mitkDICOMSegmentationIO.cpp +++ b/Modules/DICOMQI/autoload/IO/mitkDICOMSegmentationIO.cpp @@ -1,692 +1,694 @@ /*=================================================================== 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 __mitkDICOMSegmentationIO__cpp #define __mitkDICOMSegmentationIO__cpp #include "mitkDICOMSegmentationIO.h" #include "mitkDICOMQIIOMimeTypes.h" #include "mitkDICOMSegmentationConstants.h" #include #include #include #include #include #include #include #include // itk #include // dcmqi #include // us #include #include namespace mitk { DICOMSegmentationIO::DICOMSegmentationIO() : AbstractFileIO(LabelSetImage::GetStaticNameOfClass(), mitk::MitkDICOMQIIOMimeTypes::DICOMQI_MIMETYPE_NAME(), "DICOM Segmentation") { AbstractFileWriter::SetRanking(10); AbstractFileReader::SetRanking(10); this->RegisterService(); this->AddDICOMTagsToService(); } void DICOMSegmentationIO::AddDICOMTagsToService() { IDICOMTagsOfInterest *toiService = GetDicomTagsOfInterestService(); if (toiService != nullptr) { toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_SEQUENCE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_NUMBER_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_LABEL_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_ALGORITHM_TYPE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::ANATOMIC_REGION_SEQUENCE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::ANATOMIC_REGION_CODE_VALUE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::ANATOMIC_REGION_CODE_SCHEME_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::ANATOMIC_REGION_CODE_MEANING_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENTED_PROPERTY_CATEGORY_SEQUENCE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_VALUE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_SCHEME_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_MEANING_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENTED_PROPERTY_TYPE_SEQUENCE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_VALUE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_SCHEME_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_MEANING_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENTED_PROPERTY_MODIFIER_SEQUENCE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_VALUE_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_SCHEME_PATH()); toiService->AddTagOfInterest(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_MEANING_PATH()); } } IFileIO::ConfidenceLevel DICOMSegmentationIO::GetWriterConfidenceLevel() const { if (AbstractFileIO::GetWriterConfidenceLevel() == Unsupported) return Unsupported; const LabelSetImage *input = static_cast(this->GetInput()); if (input) return Supported; else return Unsupported; } void DICOMSegmentationIO::Write() { ValidateOutputLocation(); mitk::LocaleSwitch localeSwitch("C"); LocalFile localFile(this); const std::string path = localFile.GetFileName(); auto input = dynamic_cast(this->GetInput()); if (input == nullptr) mitkThrow() << "Cannot write non-image data"; // Get DICOM information from referenced image vector> dcmDatasetsSourceImage; std::unique_ptr readFileFormat(new DcmFileFormat()); try { // TODO: Generate dcmdataset witk DICOM tags from property list; ATM the source are the filepaths from the // property list mitk::StringLookupTableProperty::Pointer filesProp = dynamic_cast(input->GetProperty("referenceFiles").GetPointer()); if (filesProp.IsNull()) { mitkThrow() << "No property with dicom file path."; return; } StringLookupTable filesLut = filesProp->GetValue(); const StringLookupTable::LookupTableType &lookUpTableMap = filesLut.GetLookupTable(); for (auto it : lookUpTableMap) { const char *fileName = (it.second).c_str(); if (readFileFormat->loadFile(fileName, EXS_Unknown).good()) { std::unique_ptr readDCMDataset(readFileFormat->getAndRemoveDataset()); dcmDatasetsSourceImage.push_back(std::move(readDCMDataset)); } } } catch (const std::exception &e) { MITK_ERROR << "An error occurred while getting the dicom informations: " << e.what() << endl; return; } // Iterate over all layers. For each a dcm file will be generated for (unsigned int layer = 0; layer < input->GetNumberOfLayers(); ++layer) { vector segmentations; try { // Hack: Remove the const attribute to switch between the layer images. Normally you could get the different // layer images by input->GetLayerImage(layer) mitk::LabelSetImage *mitkLayerImage = const_cast(input); mitkLayerImage->SetActiveLayer(layer); // Cast mitk layer image to itk ImageToItk::Pointer imageToItkFilter = ImageToItk::New(); imageToItkFilter->SetInput(mitkLayerImage); // Cast from original itk type to dcmqi input itk image type typedef itk::CastImageFilter castItkImageFilterType; castItkImageFilterType::Pointer castFilter = castItkImageFilterType::New(); castFilter->SetInput(imageToItkFilter->GetOutput()); castFilter->Update(); itkInternalImageType::Pointer itkLabelImage = castFilter->GetOutput(); itkLabelImage->DisconnectPipeline(); // Iterate over all labels. For each label a segmentation image will be created const LabelSet *labelSet = input->GetLabelSet(layer); auto labelIter = labelSet->IteratorConstBegin(); // Ignore background label ++labelIter; for (; labelIter != labelSet->IteratorConstEnd(); ++labelIter) { // Thresold over the image with the given label value itk::ThresholdImageFilter::Pointer thresholdFilter = itk::ThresholdImageFilter::New(); thresholdFilter->SetInput(itkLabelImage); thresholdFilter->ThresholdOutside(labelIter->first, labelIter->first); thresholdFilter->SetOutsideValue(0); thresholdFilter->Update(); itkInternalImageType::Pointer segmentImage = thresholdFilter->GetOutput(); segmentImage->DisconnectPipeline(); segmentations.push_back(segmentImage); } } catch (const itk::ExceptionObject &e) { MITK_ERROR << e.GetDescription() << endl; return; } // Create segmentation meta information const std::string tmpMetaInfoFile = this->CreateMetaDataJsonFile(layer); MITK_INFO << "Writing image: " << path << std::endl; try { //TODO is there a better way? Interface expects a vector of raw pointer. vector rawVecDataset; for (const auto& dcmDataSet : dcmDatasetsSourceImage) rawVecDataset.push_back(dcmDataSet.get()); // Convert itk segmentation images to dicom image std::unique_ptr converter = std::make_unique(); std::unique_ptr result(converter->itkimage2dcmSegmentation(rawVecDataset, segmentations, tmpMetaInfoFile)); // Write dicom file DcmFileFormat dcmFileFormat(result.get()); std::string filePath = path.substr(0, path.find_last_of(".")); // If there is more than one layer, we have to write more than 1 dicom file if (input->GetNumberOfLayers() != 1) filePath = filePath + std::to_string(layer) + ".dcm"; else filePath = filePath + ".dcm"; dcmFileFormat.saveFile(filePath.c_str(), EXS_LittleEndianExplicit); } catch (const std::exception &e) { MITK_ERROR << "An error occurred during writing the DICOM Seg: " << e.what() << endl; return; } } // Write a dcm file for the next layer } IFileIO::ConfidenceLevel DICOMSegmentationIO::GetReaderConfidenceLevel() const { if (AbstractFileIO::GetReaderConfidenceLevel() == Unsupported) return Unsupported; const std::string fileName = this->GetLocalFileName(); DcmFileFormat dcmFileFormat; OFCondition status = dcmFileFormat.loadFile(fileName.c_str()); if (status.bad()) return Unsupported; OFString modality; if (dcmFileFormat.getDataset()->findAndGetOFString(DCM_Modality, modality).good()) { if (modality.compare("SEG") == 0) return Supported; else return Unsupported; } return Unsupported; } std::vector DICOMSegmentationIO::Read() { mitk::LocaleSwitch localeSwitch("C"); LabelSetImage::Pointer labelSetImage; std::vector result; const std::string path = this->GetLocalFileName(); MITK_INFO << "loading " << path << std::endl; if (path.empty()) mitkThrow() << "Empty filename in mitk::ItkImageIO "; try { // Get the dcm data set from file path DcmFileFormat dcmFileFormat; OFCondition status = dcmFileFormat.loadFile(path.c_str()); if (status.bad()) mitkThrow() << "Can't read the input file!"; DcmDataset *dataSet = dcmFileFormat.getDataset(); if (dataSet == nullptr) mitkThrow() << "Can't read data from input file!"; //=============================== dcmqi part ==================================== // Read the DICOM SEG images (segItkImages) and DICOM tags (metaInfo) std::unique_ptr converter = std::make_unique(); pair, string> dcmqiOutput = converter->dcmSegmentation2itkimage(dataSet); map segItkImages = dcmqiOutput.first; dcmqi::JSONSegmentationMetaInformationHandler metaInfo(dcmqiOutput.second.c_str()); metaInfo.read(); MITK_INFO << "Input " << metaInfo.getJSONOutputAsString(); //=============================================================================== // Get the label information from segment attributes for each itk image vector>::const_iterator segmentIter = metaInfo.segmentsAttributesMappingList.begin(); // For each itk image add a layer to the LabelSetImage output for (auto &element : segItkImages) { // Get the labeled image and cast it to mitkImage typedef itk::CastImageFilter castItkImageFilterType; castItkImageFilterType::Pointer castFilter = castItkImageFilterType::New(); castFilter->SetInput(element.second); castFilter->Update(); Image::Pointer layerImage; CastToMitkImage(castFilter->GetOutput(), layerImage); // Get pixel value of the label itkInternalImageType::ValueType segValue = 1; typedef itk::ImageRegionIterator IteratorType; // Iterate over the image to find the pixel value of the label IteratorType iter(element.second, element.second->GetLargestPossibleRegion()); iter.GoToBegin(); while (!iter.IsAtEnd()) { itkInputImageType::PixelType value = iter.Get(); if (value != 0) { segValue = value; break; } ++iter; } // Get Segment information map map segmentMap = (*segmentIter); map::const_iterator segmentMapIter = (*segmentIter).begin(); dcmqi::SegmentAttributes *segmentAttribute = (*segmentMapIter).second; OFString labelName; if (segmentAttribute->getSegmentedPropertyTypeCodeSequence() != nullptr) { segmentAttribute->getSegmentedPropertyTypeCodeSequence()->getCodeMeaning(labelName); if (segmentAttribute->getSegmentedPropertyTypeModifierCodeSequence() != nullptr) { OFString modifier; segmentAttribute->getSegmentedPropertyTypeModifierCodeSequence()->getCodeMeaning(modifier); labelName.append(" (").append(modifier).append(")"); } } else { labelName = std::to_string(segmentAttribute->getLabelID()).c_str(); if (labelName.empty()) labelName = "Unnamed"; } float tmp[3] = { 0.0, 0.0, 0.0 }; if (segmentAttribute->getRecommendedDisplayRGBValue() != nullptr) { tmp[0] = segmentAttribute->getRecommendedDisplayRGBValue()[0] / 255.0; tmp[1] = segmentAttribute->getRecommendedDisplayRGBValue()[1] / 255.0; tmp[2] = segmentAttribute->getRecommendedDisplayRGBValue()[2] / 255.0; } Label *newLabel = nullptr; // If labelSetImage do not exists (first image) if (labelSetImage.IsNull()) { // Initialize the labelSetImage with the read image labelSetImage = LabelSetImage::New(); labelSetImage->InitializeByLabeledImage(layerImage); // Already a label was generated, so set the information to this newLabel = labelSetImage->GetActiveLabel(labelSetImage->GetActiveLayer()); newLabel->SetName(labelName.c_str()); newLabel->SetColor(Color(tmp)); newLabel->SetValue(segValue); } else { // Add a new layer to the labelSetImage. Background label is set automatically labelSetImage->AddLayer(layerImage); // Add new label newLabel = new Label; newLabel->SetName(labelName.c_str()); newLabel->SetColor(Color(tmp)); newLabel->SetValue(segValue); labelSetImage->GetLabelSet(labelSetImage->GetActiveLayer())->AddLabel(newLabel); } // Add some more label properties this->SetLabelProperties(newLabel, segmentAttribute); ++segmentIter; } // Add some general DICOM Segmentation properties mitk::IDICOMTagsOfInterest *toiSrv = GetDicomTagsOfInterestService(); auto tagsOfInterest = toiSrv->GetTagsOfInterest(); DICOMTagPathList tagsOfInterestList; for (const auto &tag : tagsOfInterest) { tagsOfInterestList.push_back(tag.first); } mitk::DICOMDCMTKTagScanner::Pointer scanner = mitk::DICOMDCMTKTagScanner::New(); scanner->SetInputFiles({ GetInputLocation() }); scanner->AddTagPaths(tagsOfInterestList); scanner->Scan(); mitk::DICOMDatasetAccessingImageFrameList frames = scanner->GetFrameInfoList(); if (frames.empty()) { MITK_ERROR << "Error reading the DICOM Seg file" << std::endl; return result; } auto findings = ExtractPathsOfInterest(tagsOfInterestList, frames); SetProperties(labelSetImage, findings); // Set active layer to the first layer of the labelset image if (labelSetImage->GetNumberOfLayers() > 1 && labelSetImage->GetActiveLayer() != 0) labelSetImage->SetActiveLayer(0); } catch (const std::exception &e) { MITK_ERROR << "An error occurred while reading the DICOM Seg file: " << e.what(); return result; } catch (...) { MITK_ERROR << "An error occurred in dcmqi while reading the DICOM Seg file"; return result; } result.push_back(labelSetImage.GetPointer()); return result; } const std::string mitk::DICOMSegmentationIO::CreateMetaDataJsonFile(int layer) { const mitk::LabelSetImage *image = dynamic_cast(this->GetInput()); const std::string output; dcmqi::JSONSegmentationMetaInformationHandler handler; // 1. Metadata attributes that will be listed in the resulting DICOM SEG object std::string contentCreatorName; if (!image->GetPropertyList()->GetStringProperty(GeneratePropertyNameForDICOMTag(0x0070, 0x0084).c_str(), contentCreatorName)) contentCreatorName = "MITK"; handler.setContentCreatorName(contentCreatorName); std::string clinicalTrailSeriesId; if (!image->GetPropertyList()->GetStringProperty(GeneratePropertyNameForDICOMTag(0x0012, 0x0071).c_str(), clinicalTrailSeriesId)) clinicalTrailSeriesId = "Session 1"; handler.setClinicalTrialSeriesID(clinicalTrailSeriesId); std::string clinicalTrialTimePointID; if (!image->GetPropertyList()->GetStringProperty(GeneratePropertyNameForDICOMTag(0x0012, 0x0050).c_str(), clinicalTrialTimePointID)) clinicalTrialTimePointID = "0"; handler.setClinicalTrialTimePointID(clinicalTrialTimePointID); std::string clinicalTrialCoordinatingCenterName = ""; if (!image->GetPropertyList()->GetStringProperty(GeneratePropertyNameForDICOMTag(0x0012, 0x0060).c_str(), clinicalTrialCoordinatingCenterName)) clinicalTrialCoordinatingCenterName = "Unknown"; handler.setClinicalTrialCoordinatingCenterName(clinicalTrialCoordinatingCenterName); std::string seriesDescription; if (!image->GetPropertyList()->GetStringProperty("name", seriesDescription)) seriesDescription = "MITK Segmentation"; handler.setSeriesDescription(seriesDescription); handler.setSeriesNumber("0" + std::to_string(layer)); handler.setInstanceNumber("1"); handler.setBodyPartExamined(""); const LabelSet *labelSet = image->GetLabelSet(layer); auto labelIter = labelSet->IteratorConstBegin(); + MITK_INFO << labelIter->second->GetValue(); // Ignore background label ++labelIter; + MITK_INFO << labelIter->second->GetValue(); for (; labelIter != labelSet->IteratorConstEnd(); ++labelIter) { const Label *label = labelIter->second; if (label != nullptr) { TemporoSpatialStringProperty *segmentNumberProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_NUMBER_PATH()).c_str())); TemporoSpatialStringProperty *segmentLabelProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_LABEL_PATH()).c_str())); TemporoSpatialStringProperty *algorithmTypeProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_ALGORITHM_TYPE_PATH()).c_str())); TemporoSpatialStringProperty *segmentCategoryCodeValueProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_VALUE_PATH()).c_str())); TemporoSpatialStringProperty *segmentCategoryCodeSchemeProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_SCHEME_PATH()).c_str())); TemporoSpatialStringProperty *segmentCategoryCodeMeaningProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_MEANING_PATH()).c_str())); TemporoSpatialStringProperty *segmentTypeCodeValueProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_VALUE_PATH()).c_str())); TemporoSpatialStringProperty *segmentTypeCodeSchemeProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_SCHEME_PATH()).c_str())); TemporoSpatialStringProperty *segmentTypeCodeMeaningProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_MEANING_PATH()).c_str())); TemporoSpatialStringProperty *segmentModifierCodeValueProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_VALUE_PATH()).c_str())); TemporoSpatialStringProperty *segmentModifierCodeSchemeProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_SCHEME_PATH()).c_str())); TemporoSpatialStringProperty *segmentModifierCodeMeaningProp = dynamic_cast(label->GetProperty( mitk::DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_MEANING_PATH()).c_str())); dcmqi::SegmentAttributes *segmentAttribute = nullptr; if (segmentNumberProp->GetValue() == "") { MITK_ERROR << "Something went wrong with the label ID."; } else { int labelId = std::stoi(segmentNumberProp->GetValue()); segmentAttribute = handler.createAndGetNewSegment(labelId); } if (segmentAttribute != nullptr) { segmentAttribute->setSegmentDescription(segmentLabelProp->GetValueAsString()); segmentAttribute->setSegmentAlgorithmType(algorithmTypeProp->GetValueAsString()); segmentAttribute->setSegmentAlgorithmName("MITK Segmentation"); if (segmentCategoryCodeValueProp != nullptr && segmentCategoryCodeSchemeProp != nullptr && segmentCategoryCodeMeaningProp != nullptr) segmentAttribute->setSegmentedPropertyCategoryCodeSequence( segmentCategoryCodeValueProp->GetValueAsString(), segmentCategoryCodeSchemeProp->GetValueAsString(), segmentCategoryCodeMeaningProp->GetValueAsString()); else // some default values segmentAttribute->setSegmentedPropertyCategoryCodeSequence( "M-01000", "SRT", "Morphologically Altered Structure"); if (segmentTypeCodeValueProp != nullptr && segmentTypeCodeSchemeProp != nullptr && segmentTypeCodeMeaningProp != nullptr) { segmentAttribute->setSegmentedPropertyTypeCodeSequence(segmentTypeCodeValueProp->GetValueAsString(), segmentTypeCodeSchemeProp->GetValueAsString(), segmentTypeCodeMeaningProp->GetValueAsString()); handler.setBodyPartExamined(segmentTypeCodeMeaningProp->GetValueAsString()); } else { // some default values segmentAttribute->setSegmentedPropertyTypeCodeSequence("M-03000", "SRT", "Mass"); handler.setBodyPartExamined("Mass"); } if (segmentModifierCodeValueProp != nullptr && segmentModifierCodeSchemeProp != nullptr && segmentModifierCodeMeaningProp != nullptr) segmentAttribute->setSegmentedPropertyTypeModifierCodeSequence( segmentModifierCodeValueProp->GetValueAsString(), segmentModifierCodeSchemeProp->GetValueAsString(), segmentModifierCodeMeaningProp->GetValueAsString()); Color color = label->GetColor(); segmentAttribute->setRecommendedDisplayRGBValue(color[0] * 255, color[1] * 255, color[2] * 255); } } } return handler.getJSONOutputAsString(); } void mitk::DICOMSegmentationIO::SetLabelProperties(mitk::Label *label, dcmqi::SegmentAttributes *segmentAttribute) { // Segment Number:Identification number of the segment.The value of Segment Number(0062, 0004) shall be unique // within the Segmentation instance in which it is created label->SetProperty(DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_NUMBER_PATH()).c_str(), TemporoSpatialStringProperty::New(std::to_string(label->GetValue()))); // Segment Label: User-defined label identifying this segment. label->SetProperty(DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_LABEL_PATH()).c_str(), TemporoSpatialStringProperty::New(label->GetName())); // Segment Algorithm Type: Type of algorithm used to generate the segment. if (!segmentAttribute->getSegmentAlgorithmType().empty()) label->SetProperty(DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_ALGORITHM_TYPE_PATH()).c_str(), TemporoSpatialStringProperty::New(segmentAttribute->getSegmentAlgorithmType())); // Add Segmented Property Category Code Sequence tags auto categoryCodeSequence = segmentAttribute->getSegmentedPropertyCategoryCodeSequence(); if (categoryCodeSequence != nullptr) { OFString codeValue; // (0008,0100) Code Value categoryCodeSequence->getCodeValue(codeValue); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_VALUE_PATH()).c_str(), TemporoSpatialStringProperty::New(codeValue.c_str())); OFString codeScheme; // (0008,0102) Coding Scheme Designator categoryCodeSequence->getCodingSchemeDesignator(codeScheme); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_SCHEME_PATH()).c_str(), TemporoSpatialStringProperty::New(codeScheme.c_str())); OFString codeMeaning; // (0008,0104) Code Meaning categoryCodeSequence->getCodeMeaning(codeMeaning); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_MEANING_PATH()).c_str(), TemporoSpatialStringProperty::New(codeMeaning.c_str())); } // Add Segmented Property Type Code Sequence tags auto typeCodeSequence = segmentAttribute->getSegmentedPropertyTypeCodeSequence(); if (typeCodeSequence != nullptr) { OFString codeValue; // (0008,0100) Code Value typeCodeSequence->getCodeValue(codeValue); label->SetProperty(DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_VALUE_PATH()).c_str(), TemporoSpatialStringProperty::New(codeValue.c_str())); OFString codeScheme; // (0008,0102) Coding Scheme Designator typeCodeSequence->getCodingSchemeDesignator(codeScheme); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_SCHEME_PATH()).c_str(), TemporoSpatialStringProperty::New(codeScheme.c_str())); OFString codeMeaning; // (0008,0104) Code Meaning typeCodeSequence->getCodeMeaning(codeMeaning); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_MEANING_PATH()).c_str(), TemporoSpatialStringProperty::New(codeMeaning.c_str())); } // Add Segmented Property Type Modifier Code Sequence tags auto modifierCodeSequence = segmentAttribute->getSegmentedPropertyTypeModifierCodeSequence(); if (modifierCodeSequence != nullptr) { OFString codeValue; // (0008,0100) Code Value modifierCodeSequence->getCodeValue(codeValue); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_VALUE_PATH()).c_str(), TemporoSpatialStringProperty::New(codeValue.c_str())); OFString codeScheme; // (0008,0102) Coding Scheme Designator modifierCodeSequence->getCodingSchemeDesignator(codeScheme); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_SCHEME_PATH()).c_str(), TemporoSpatialStringProperty::New(codeScheme.c_str())); OFString codeMeaning; // (0008,0104) Code Meaning modifierCodeSequence->getCodeMeaning(codeMeaning); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_MODIFIER_CODE_MEANING_PATH()).c_str(), TemporoSpatialStringProperty::New(codeMeaning.c_str())); } // Add Atomic RegionSequence tags auto atomicRegionSequence = segmentAttribute->getAnatomicRegionSequence(); if (atomicRegionSequence != nullptr) { OFString codeValue; // (0008,0100) Code Value atomicRegionSequence->getCodeValue(codeValue); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::ANATOMIC_REGION_CODE_VALUE_PATH()).c_str(), TemporoSpatialStringProperty::New(codeValue.c_str())); OFString codeScheme; // (0008,0102) Coding Scheme Designator atomicRegionSequence->getCodingSchemeDesignator(codeScheme); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::ANATOMIC_REGION_CODE_SCHEME_PATH()).c_str(), TemporoSpatialStringProperty::New(codeScheme.c_str())); OFString codeMeaning; // (0008,0104) Code Meaning atomicRegionSequence->getCodeMeaning(codeMeaning); label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::ANATOMIC_REGION_CODE_MEANING_PATH()).c_str(), TemporoSpatialStringProperty::New(codeMeaning.c_str())); } } DICOMSegmentationIO *DICOMSegmentationIO::IOClone() const { return new DICOMSegmentationIO(*this); } } // namespace #endif //__mitkDICOMSegmentationIO__cpp diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp index a19f7a97e2..74a0ede0ac 100644 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp @@ -1,110 +1,72 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include - -template -mitk::OdfImage::Pointer TemplatedConvertShImage(mitk::ShImage::Pointer mitkImage) -{ - typedef itk::ShToOdfImageFilter< float, ShOrder > ShConverterType; - typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); - mitk::CastToItkImage(mitkImage, itkvol); - - typename ShConverterType::Pointer converter = ShConverterType::New(); - converter->SetInput(itkvol); - converter->Update(); - - mitk::OdfImage::Pointer image = mitk::OdfImage::New(); - image->InitializeByItk( converter->GetOutput() ); - image->SetVolume( converter->GetOutput()->GetBufferPointer() ); - return image; -} +#include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("ShToOdfImage"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Calculate discrete ODF image from SH coefficient image"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputFile, "Output:", "output image", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); try { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); mitk::ShImage::Pointer source = mitk::IOUtil::Load(imageName, &functor); - mitk::OdfImage::Pointer out_image = nullptr; - switch (source->ShOrder()) - { - case 2: - out_image = TemplatedConvertShImage<2>(source); - break; - case 4: - out_image = TemplatedConvertShImage<4>(source); - break; - case 6: - out_image = TemplatedConvertShImage<6>(source); - break; - case 8: - out_image = TemplatedConvertShImage<8>(source); - break; - case 10: - out_image = TemplatedConvertShImage<10>(source); - break; - case 12: - out_image = TemplatedConvertShImage<12>(source); - break; - } - + mitk::OdfImage::Pointer out_image = mitk::convert::GetOdfFromShImage(source); if (out_image.IsNotNull()) mitk::IOUtil::Save(out_image, outImage); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp index b5cf08084e..3737e98a02 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp @@ -1,170 +1,133 @@ /*=================================================================== 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 -typename itk::ShCoefficientImageImporter< float, shOrder >::OdfImageType::Pointer TemplatedConvertShCoeffs(mitk::Image::Pointer mitkImg) -{ - typedef itk::ShToOdfImageFilter< float, shOrder > ShConverterType; - - typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); - mitk::CastToItkImage(mitkImg, itkvol); - - typename ShConverterType::Pointer converter = ShConverterType::New(); - converter->SetInput(itkvol); - converter->Update(); - - return converter->GetOutput(); -} +#include /*! \brief Perform global fiber tractography (Gibbs tractography) */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Gibbs Tracking"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform global fiber tractography (Gibbs tractography)"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input:", "input image (tensor, ODF or SH-coefficient image)", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputFile, "Output:", "output tractogram", us::Any(), false); parser.addArgument("parameters", "", mitkCommandLineParser::InputFile, "Parameters:", "parameter file (.gtp)", us::Any(), false); parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "Mask:", "binary mask image"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string paramFileName = us::any_cast(parsedArgs["parameters"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { // instantiate gibbs tracker typedef itk::Vector OdfVectorType; typedef itk::Image ItkOdfImageType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); // load input image mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(inFileName, &functor); // try to cast to Odf image if( dynamic_cast(mitkImage.GetPointer()) ) { mitk::OdfImage::Pointer mitkOdfImage = dynamic_cast(mitkImage.GetPointer()); ItkOdfImageType::Pointer itk_odf = ItkOdfImageType::New(); mitk::CastToItkImage(mitkOdfImage, itk_odf); gibbsTracker->SetOdfImage(itk_odf.GetPointer()); } else if( dynamic_cast(mitkImage.GetPointer()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(mitkImage.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); gibbsTracker->SetTensorImage(itk_dti); } else if ( dynamic_cast(mitkImage.GetPointer()) ) { mitk::ShImage::Pointer shImage = dynamic_cast(mitkImage.GetPointer()); - - switch (shImage->ShOrder()) - { - case 2: - gibbsTracker->SetOdfImage(TemplatedConvertShCoeffs<2>(mitkImage)); - break; - case 4: - gibbsTracker->SetOdfImage(TemplatedConvertShCoeffs<4>(mitkImage)); - break; - case 6: - gibbsTracker->SetOdfImage(TemplatedConvertShCoeffs<6>(mitkImage)); - break; - case 8: - gibbsTracker->SetOdfImage(TemplatedConvertShCoeffs<8>(mitkImage)); - break; - case 10: - gibbsTracker->SetOdfImage(TemplatedConvertShCoeffs<10>(mitkImage)); - break; - case 12: - gibbsTracker->SetOdfImage(TemplatedConvertShCoeffs<12>(mitkImage)); - break; - default: - std::cout << "SH-order " << shImage->ShOrder() << " not supported"; - } + gibbsTracker->SetOdfImage(mitk::convert::GetItkOdfFromShImage(shImage)); } else return EXIT_FAILURE; // global tracking if (parsedArgs.count("mask")) { typedef itk::Image MaskImgType; mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::Load(us::any_cast(parsedArgs["mask"])); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); gibbsTracker->SetMaskImage(itk_mask); } gibbsTracker->SetDuplicateImage(false); gibbsTracker->SetLoadParameterFile( paramFileName ); // gibbsTracker->SetLutPath( "" ); gibbsTracker->Update(); mitk::FiberBundle::Pointer mitkFiberBundle = mitk::FiberBundle::New(gibbsTracker->GetFiberBundle()); mitkFiberBundle->SetReferenceGeometry(mitkImage->GetGeometry()); mitk::IOUtil::Save(mitkFiberBundle, outFileName ); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp index e84ca3ed06..a885a256ec 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp @@ -1,573 +1,561 @@ -/*=================================================================== +/*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle/probability map", us::Any(), false); parser.addArgument("algorithm", "", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF)", us::Any(), false); parser.endGroup(); parser.beginGroup("2. Seeding:"); parser.addArgument("seeds", "", mitkCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_image", "", mitkCommandLineParser::String, "Seed image:", "mask image defining seed voxels", us::Any()); parser.addArgument("trials_per_seed", "", mitkCommandLineParser::Int, "Max. trials per seed:", "try each seed N times until a valid streamline is obtained (only for probabilistic tractography)", 10); parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded", -1); parser.endGroup(); parser.beginGroup("3. Tractography constraints:"); parser.addArgument("tracking_mask", "", mitkCommandLineParser::String, "Mask image:", "streamlines leaving the mask will stop immediately", us::Any()); parser.addArgument("stop_image", "", mitkCommandLineParser::String, "Stop ROI image:", "streamlines entering the mask will stop immediately", us::Any()); parser.addArgument("exclusion_image", "", mitkCommandLineParser::String, "Exclusion ROI image:", "streamlines entering the mask will be discarded", us::Any()); parser.addArgument("ep_constraint", "", mitkCommandLineParser::String, "Endpoint constraint:", "determines which fibers are accepted based on their endpoint location - options are NONE, EPS_IN_TARGET, EPS_IN_TARGET_LABELDIFF, EPS_IN_SEED_AND_TARGET, MIN_ONE_EP_IN_TARGET, ONE_EP_IN_TARGET and NO_EP_IN_TARGET", us::Any()); parser.addArgument("target_image", "", mitkCommandLineParser::String, "Target ROI image:", "effact depends on the chosen endpoint constraint (option ep_constraint)", us::Any()); parser.endGroup(); parser.beginGroup("4. Streamline integration parameters:"); parser.addArgument("sharpen_odfs", "", mitkCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkCommandLineParser::Float, "ODF Cutoff:", "threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.0); parser.addArgument("step_size", "", mitkCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("min_tract_length", "", mitkCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("angular_threshold", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size, minimum 15°)"); parser.addArgument("loop_check", "", mitkCommandLineParser::Float, "Check for loops:", "threshold on angular stdev over the last 4 voxel lengths"); parser.endGroup(); parser.beginGroup("5. Tractography prior:"); parser.addArgument("prior_image", "", mitkCommandLineParser::String, "Peak prior:", "tractography prior in thr for of a peak image", us::Any()); parser.addArgument("prior_weight", "", mitkCommandLineParser::Float, "Prior weight", "weighting factor between prior and data.", 0.5); parser.addArgument("restrict_to_prior", "", mitkCommandLineParser::Bool, "Restrict to prior:", "restrict tractography to regions where the prior is valid."); parser.addArgument("new_directions_from_prior", "", mitkCommandLineParser::Bool, "New directios from prior:", "the prior can create directions where there are none in the data."); parser.endGroup(); parser.beginGroup("6. Neighborhood sampling:"); parser.addArgument("num_samples", "", mitkCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.endGroup(); parser.beginGroup("7. Tensor tractography specific:"); parser.addArgument("tend_f", "", mitkCommandLineParser::Float, "Weight f", "weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkCommandLineParser::Float, "Weight g", "weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); parser.endGroup(); parser.beginGroup("8. Random forest tractography specific:"); parser.addArgument("forest", "", mitkCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any()); parser.addArgument("use_sh_features", "", mitkCommandLineParser::Bool, "Use SH features:", "use SH features"); parser.endGroup(); parser.beginGroup("9. Additional input:"); parser.addArgument("additional_images", "", mitkCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional features for RF tractography)", us::Any()); parser.endGroup(); parser.beginGroup("10. Misc:"); parser.addArgument("flip_x", "", mitkCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); parser.addArgument("no_data_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate input data:", "don't interpolate input image values"); parser.addArgument("no_mask_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate masks:", "don't interpolate mask image values"); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "compress output fibers using the given error threshold (in mm)"); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["i"]); std::string outFile = us::any_cast(parsedArgs["o"]); std::string algorithm = us::any_cast(parsedArgs["algorithm"]); std::string prior_image = ""; if (parsedArgs.count("prior_image")) prior_image = us::any_cast(parsedArgs["prior_image"]); float prior_weight = 0.5; if (parsedArgs.count("prior_weight")) prior_weight = us::any_cast(parsedArgs["prior_weight"]); bool restrict_to_prior = false; if (parsedArgs.count("restrict_to_prior")) restrict_to_prior = us::any_cast(parsedArgs["restrict_to_prior"]); bool new_directions_from_prior = false; if (parsedArgs.count("new_directions_from_prior")) new_directions_from_prior = us::any_cast(parsedArgs["new_directions_from_prior"]); bool sharpen_odfs = false; if (parsedArgs.count("sharpen_odfs")) sharpen_odfs = us::any_cast(parsedArgs["sharpen_odfs"]); bool interpolate = true; if (parsedArgs.count("no_data_interpolation")) interpolate = !us::any_cast(parsedArgs["no_data_interpolation"]); bool mask_interpolation = true; if (parsedArgs.count("no_mask_interpolation")) interpolate = !us::any_cast(parsedArgs["no_mask_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); bool use_stop_votes = false; if (parsedArgs.count("use_stop_votes")) use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); bool use_only_forward_samples = false; if (parsedArgs.count("use_only_forward_samples")) use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); bool flip_x = false; if (parsedArgs.count("flip_x")) flip_x = us::any_cast(parsedArgs["flip_x"]); bool flip_y = false; if (parsedArgs.count("flip_y")) flip_y = us::any_cast(parsedArgs["flip_y"]); bool flip_z = false; if (parsedArgs.count("flip_z")) flip_z = us::any_cast(parsedArgs["flip_z"]); bool apply_image_rotation = false; if (parsedArgs.count("apply_image_rotation")) apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); float min_tract_length = 20; if (parsedArgs.count("min_tract_length")) min_tract_length = us::any_cast(parsedArgs["min_tract_length"]); float loop_check = -1; if (parsedArgs.count("loop_check")) loop_check = us::any_cast(parsedArgs["loop_check"]); std::string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); std::string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); std::string seedFile = ""; if (parsedArgs.count("seed_image")) seedFile = us::any_cast(parsedArgs["seed_image"]); std::string targetFile = ""; if (parsedArgs.count("target_image")) targetFile = us::any_cast(parsedArgs["target_image"]); std::string exclusionFile = ""; if (parsedArgs.count("exclusion_image")) exclusionFile = us::any_cast(parsedArgs["exclusion_image"]); std::string stopFile = ""; if (parsedArgs.count("stop_image")) stopFile = us::any_cast(parsedArgs["stop_image"]); std::string ep_constraint = "NONE"; if (parsedArgs.count("ep_constraint")) ep_constraint = us::any_cast(parsedArgs["ep_constraint"]); - float cutoff = 0.1; + float cutoff = 0.1f; if (parsedArgs.count("cutoff")) cutoff = us::any_cast(parsedArgs["cutoff"]); float odf_cutoff = 0.0; if (parsedArgs.count("odf_cutoff")) odf_cutoff = us::any_cast(parsedArgs["odf_cutoff"]); float stepsize = -1; if (parsedArgs.count("step_size")) stepsize = us::any_cast(parsedArgs["step_size"]); float sampling_distance = -1; if (parsedArgs.count("sampling_distance")) sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); int num_samples = 0; if (parsedArgs.count("num_samples")) num_samples = us::any_cast(parsedArgs["num_samples"]); int num_seeds = 1; if (parsedArgs.count("seeds")) num_seeds = us::any_cast(parsedArgs["seeds"]); unsigned int trials_per_seed = 10; if (parsedArgs.count("trials_per_seed")) - trials_per_seed = us::any_cast(parsedArgs["trials_per_seed"]); + trials_per_seed = us::any_cast(parsedArgs["trials_per_seed"]); float tend_f = 1; if (parsedArgs.count("tend_f")) tend_f = us::any_cast(parsedArgs["tend_f"]); float tend_g = 0; if (parsedArgs.count("tend_g")) tend_g = us::any_cast(parsedArgs["tend_g"]); float angular_threshold = -1; if (parsedArgs.count("angular_threshold")) angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); - unsigned int max_tracts = -1; + int max_tracts = -1; if (parsedArgs.count("max_tracts")) max_tracts = us::any_cast(parsedArgs["max_tracts"]); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (ext != ".fib" && ext != ".trk") { MITK_INFO << "Output file format not supported. Use one of .fib, .trk, .nii, .nii.gz, .nrrd"; return EXIT_FAILURE; } // LOAD DATASETS mitkCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkFloatImgType; - MITK_INFO << "loading input"; - std::vector< mitk::Image::Pointer > input_images; - for (unsigned int i=0; i(input_files.at(i)); - input_images.push_back(mitkImage); - } - ItkFloatImgType::Pointer mask = nullptr; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = mitk::IOUtil::Load(maskFile); mask = ItkFloatImgType::New(); mitk::CastToItkImage(img, mask); } ItkFloatImgType::Pointer seed = nullptr; if (!seedFile.empty()) { MITK_INFO << "loading seed ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(seedFile); seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); } ItkFloatImgType::Pointer stop = nullptr; if (!stopFile.empty()) { MITK_INFO << "loading stop ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(stopFile); stop = ItkFloatImgType::New(); mitk::CastToItkImage(img, stop); } ItkFloatImgType::Pointer target = nullptr; if (!targetFile.empty()) { MITK_INFO << "loading target ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(targetFile); target = ItkFloatImgType::New(); mitk::CastToItkImage(img, target); } ItkFloatImgType::Pointer exclusion = nullptr; if (!exclusionFile.empty()) { MITK_INFO << "loading exclusion ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(exclusionFile); exclusion = ItkFloatImgType::New(); mitk::CastToItkImage(img, exclusion); } MITK_INFO << "loading additional images"; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(file); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// -// omp_set_num_threads(1); - - if (algorithm == "ProbTensor") - { - typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(input_images.at(0)); - caster->Update(); - mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkTensorImg = caster->GetOutput(); - - typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; - FilterType::Pointer filter = FilterType::New(); - filter->SetInput( itkTensorImg ); - filter->Update(); - - mitk::Image::Pointer image = mitk::Image::New(); - FilterType::OutputImageType::Pointer outimg = filter->GetOutput(); - image->InitializeByItk( outimg.GetPointer() ); - image->SetVolume( outimg->GetBufferPointer() ); - - input_images.clear(); - input_images.push_back(image); - - sharpen_odfs = true; - odf_cutoff = 0; - } + // omp_set_num_threads(1); typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if (!prior_image.empty()) { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); mitk::PeakImage::Pointer priorImage = mitk::IOUtil::Load(prior_image, &functor); if (priorImage.IsNull()) { MITK_INFO << "Only peak images are supported as prior at the moment!"; return EXIT_FAILURE; } mitk::TrackingDataHandler* priorhandler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(priorImage); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); dynamic_cast(priorhandler)->SetPeakImage(itkImg); dynamic_cast(priorhandler)->SetPeakThreshold(0.0); dynamic_cast(priorhandler)->SetInterpolate(interpolate); dynamic_cast(priorhandler)->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); tracker->SetTrackingPriorHandler(priorhandler); tracker->SetTrackingPriorWeight(prior_weight); tracker->SetTrackingPriorAsMask(restrict_to_prior); tracker->SetIntroduceDirectionsFromPrior(new_directions_from_prior); } mitk::TrackingDataHandler* handler; if (algorithm == "DetRF" || algorithm == "ProbRF") { mitk::TractographyForest::Pointer forest = mitk::IOUtil::Load(forestFile); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + auto input = mitk::IOUtil::Load(input_files.at(0), &functor); + if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->SetForest(forest); - dynamic_cast*>(handler)->AddDwi(input_images.at(0)); + dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->SetForest(forest); - dynamic_cast*>(handler)->AddDwi(input_images.at(0)); + dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } if (algorithm == "ProbRF") handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); } else if (algorithm == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); - typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(input_images.at(0)); - caster->Update(); - mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); + MITK_INFO << "loading input peak image"; + mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(input_files.at(0)); + mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = mitk::convert::GetItkPeakFromPeakImage(mitkImage); dynamic_cast(handler)->SetPeakImage(itkImg); dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); dynamic_cast(handler)->SetPeakThreshold(cutoff); } else if (algorithm == "DetTensor") { handler = new mitk::TrackingHandlerTensor(); - for (auto input_image : input_images) + MITK_INFO << "loading input tensor images"; + std::vector< mitk::Image::Pointer > input_images; + for (unsigned int i=0; i CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(input_image); - caster->Update(); - mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); - dynamic_cast(handler)->AddTensorImage(itkImg); + mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(input_files.at(i)); + mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::convert::GetItkTensorFromTensorImage(mitkImage); + dynamic_cast(handler)->AddTensorImage(itkImg.GetPointer()); } dynamic_cast(handler)->SetFaThreshold(cutoff); dynamic_cast(handler)->SetF(tend_f); dynamic_cast(handler)->SetG(tend_g); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (algorithm == "DetODF" || algorithm == "ProbODF" || algorithm == "ProbTensor") { handler = new mitk::TrackingHandlerOdf(); - typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(input_images.at(0)); - caster->Update(); - mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); - dynamic_cast(handler)->SetOdfImage(itkImg); + mitk::OdfImage::ItkOdfImageType::Pointer itkImg = nullptr; + + if (algorithm == "ProbTensor") + { + MITK_INFO << "Converting Tensor to ODF image"; + auto input = mitk::IOUtil::Load(input_files.at(0)); + itkImg = mitk::convert::GetItkOdfFromTensorImage(input); + sharpen_odfs = true; + odf_cutoff = 0; + } + else + { + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image", "ODF Image"}, {}); + auto input = mitk::IOUtil::Load(input_files.at(0), &functor)[0]; + if (dynamic_cast(input.GetPointer())) + { + MITK_INFO << "Converting SH to ODF image"; + mitk::ShImage::Pointer mitkImg = dynamic_cast(input.GetPointer()); + itkImg = mitk::convert::GetItkOdfFromShImage(mitkImg); + } + else if (dynamic_cast(input.GetPointer())) + { + mitk::OdfImage::Pointer mitkImg = dynamic_cast(input.GetPointer()); + itkImg = mitk::convert::GetItkOdfFromOdfImage(mitkImg); + } + } + dynamic_cast(handler)->SetOdfImage(itkImg); dynamic_cast(handler)->SetGfaThreshold(cutoff); dynamic_cast(handler)->SetOdfThreshold(odf_cutoff); dynamic_cast(handler)->SetSharpenOdfs(sharpen_odfs); if (algorithm == "ProbODF" || algorithm == "ProbTensor") dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); if (algorithm == "ProbTensor") dynamic_cast(handler)->SetIsOdfFromTensor(true); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { MITK_INFO << "Unknown tractography algorithm (" + algorithm+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } handler->SetInterpolate(interpolate); handler->SetFlipX(flip_x); handler->SetFlipY(flip_y); handler->SetFlipZ(flip_z); if (ep_constraint=="NONE") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); else if (ep_constraint=="EPS_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); else if (ep_constraint=="EPS_IN_TARGET_LABELDIFF") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); else if (ep_constraint=="EPS_IN_SEED_AND_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); else if (ep_constraint=="MIN_ONE_EP_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); else if (ep_constraint=="ONE_EP_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); else if (ep_constraint=="NO_EP_IN_TARGET") tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); MITK_INFO << "Tractography algorithm: " << algorithm; tracker->SetInterpolateMasks(mask_interpolation); tracker->SetNumberOfSamples(num_samples); tracker->SetAngularThreshold(angular_threshold); tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetTargetRegions(target); tracker->SetExclusionRegions(exclusion); tracker->SetSeedsPerVoxel(num_seeds); tracker->SetStepSize(stepsize); tracker->SetSamplingDistance(sampling_distance); tracker->SetUseStopVotes(use_stop_votes); tracker->SetOnlyForwardSamples(use_only_forward_samples); tracker->SetLoopCheck(loop_check); tracker->SetMaxNumTracts(max_tracts); tracker->SetTrialsPerSeed(trials_per_seed); tracker->SetTrackingHandler(handler); if (ext != ".fib" && ext != ".trk") tracker->SetUseOutputProbabilityMap(true); tracker->SetMinTractLength(min_tract_length); tracker->Update(); if (ext == ".fib" || ext == ".trk") { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/mitkShImage.h b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/mitkShImage.h index a8a73aaf54..9aed1448e1 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 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(); + unsigned int ShOrder(); + unsigned int NumCoefficients(); protected: ShImage(); ~ShImage() override; template void Construct() const; - void PrintSelf(std::ostream &os, itk::Indent indent) const override; + void PrintSelf(std::ostream &os, itk::Indent indent) const override; mutable mitk::Image::Pointer m_RgbImage; - int m_ShOrder; - int m_NumCoefficients; + unsigned int m_ShOrder; + unsigned 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 a4862cbc3f..07d93279d8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/mitkDiffusionFunctionCollection.h @@ -1,148 +1,164 @@ /*=================================================================== 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 #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] = static_cast(point[0]); itkPoint[1] = static_cast(point[1]); itkPoint[2] = static_cast(point[2]); return itkPoint; } template< class TType > static itk::Point GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = static_cast(point[0]); itkPoint[1] = static_cast(point[1]); itkPoint[2] = static_cast(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; + float 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 convert +{ +public: + static mitk::OdfImage::ItkOdfImageType::Pointer GetItkOdfFromTensorImage(mitk::Image::Pointer mitkImage); + static mitk::OdfImage::Pointer GetOdfFromTensorImage(mitk::Image::Pointer mitkImage); + static mitk::OdfImage::ItkOdfImageType::Pointer GetItkOdfFromShImage(mitk::ShImage::Pointer mitkImage); + static mitk::OdfImage::Pointer GetOdfFromShImage(mitk::ShImage::Pointer mitkImage); + static mitk::OdfImage::ItkOdfImageType::Pointer GetItkOdfFromOdfImage(mitk::OdfImage::Pointer mitkImage); + static mitk::TensorImage::ItkTensorImageType::Pointer GetItkTensorFromTensorImage(mitk::Image::Pointer mitkImage); + static mitk::PeakImage::ItkPeakImageType::Pointer GetItkPeakFromPeakImage(mitk::Image::Pointer mitkImage); +}; + 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/IODataStructures/mitkShImage.cpp b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/mitkShImage.cpp index cb68f88425..90d273e9a8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/mitkShImage.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/mitkShImage.cpp @@ -1,156 +1,144 @@ /*=================================================================== 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 "mitkShImage.h" #include "mitkImageCast.h" #include "itkImage.h" #include "mitkImageVtkAccessor.h" #include #include "itkShToRgbImageFilter.h" mitk::ShImage::ShImage() : Image() , m_ShOrder(0) , m_NumCoefficients(0) { m_RgbImage = nullptr; // not needed anymore as soon as all diffusion images are identified via properties anyway this->SetProperty("IsShImage", mitk::BoolProperty::New(true)); } mitk::ShImage::~ShImage() { } -int mitk::ShImage::NumCoefficients() +unsigned int mitk::ShImage::NumCoefficients() { if (m_NumCoefficients==0 || m_ShOrder==0) { - m_NumCoefficients = this->m_ImageDescriptor->GetChannelTypeById(0).GetNumberOfComponents(); - int c=3, d=2-2*m_NumCoefficients; - double D = c*c-4*d; + m_NumCoefficients = static_cast(this->m_ImageDescriptor->GetChannelTypeById(0).GetNumberOfComponents()); + int c=3, d=2-2*static_cast(m_NumCoefficients); + int D = c*c-4*d; if (D>0) { - m_ShOrder = (-c+sqrt(D))/2.0; - if (m_ShOrder<0) - m_ShOrder = (-c-sqrt(D))/2.0; + int s = (-c+static_cast(sqrt(D)))/2; + if (s<0) + s = (-c-static_cast(sqrt(D)))/2; + m_ShOrder = static_cast(s); } else if (D==0) - m_ShOrder = -c/2.0; + m_ShOrder = static_cast(-c/2); } return m_NumCoefficients; } -int mitk::ShImage::ShOrder() +unsigned int mitk::ShImage::ShOrder() { - if (m_NumCoefficients==0 || m_ShOrder==0) - { - m_NumCoefficients = this->m_ImageDescriptor->GetChannelTypeById(0).GetNumberOfComponents(); - int c=3, d=2-2*m_NumCoefficients; - double D = c*c-4*d; - if (D>0) - { - m_ShOrder = (-c+sqrt(D))/2.0; - if (m_ShOrder<0) - m_ShOrder = (-c-sqrt(D))/2.0; - } - else if (D==0) - m_ShOrder = -c/2.0; - } + NumCoefficients(); return m_ShOrder; } vtkImageData* mitk::ShImage::GetVtkImageData(int t, int n) { if(m_RgbImage.IsNull()) { ShOrder(); NumCoefficients(); ConstructRgbImage(); } return m_RgbImage->GetVtkImageData(t,n); } const vtkImageData*mitk::ShImage::GetVtkImageData(int t, int n) const { if(m_RgbImage.IsNull()) { ConstructRgbImage(); } return m_RgbImage->GetVtkImageData(t,n); } template void mitk::ShImage::Construct() const { typedef itk::Image,3> ImageType; typedef itk::ShToRgbImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); typename ImageType::Pointer itkvol = ImageType::New(); mitk::CastToItkImage(this, itkvol); filter->SetInput(itkvol); filter->Update(); itk::Image,3>::Pointer tmp = filter->GetOutput(); m_RgbImage = mitk::Image::New(); m_RgbImage->InitializeByItk( tmp.GetPointer() ); m_RgbImage->SetVolume( tmp->GetBufferPointer() ); } void mitk::ShImage::ConstructRgbImage() const { switch (this->m_ImageDescriptor->GetChannelTypeById(0).GetNumberOfComponents()) { case 6: Construct<2>(); break; case 15: Construct<4>(); break; case 28: Construct<6>(); break; case 45: Construct<8>(); break; case 66: Construct<10>(); break; case 91: Construct<12>(); break; default : mitkThrow() << "SH order larger 12 not supported"; } } const vtkImageData* mitk::ShImage::GetNonRgbVtkImageData(int t, int n) const { return Superclass::GetVtkImageData(t,n); } vtkImageData* mitk::ShImage::GetNonRgbVtkImageData(int t, int n) { return Superclass::GetVtkImageData(t,n); } void mitk::ShImage::PrintSelf(std::ostream &os, itk::Indent indent) const { os << indent << "Spherical harmonics order: " << m_ShOrder << std::endl; os << indent << "Number of coefficients: " << m_NumCoefficients << std::endl; Superclass::PrintSelf(os, indent); } diff --git a/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp index 2279e3df08..d6174ebdb3 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/mitkDiffusionFunctionCollection.cpp @@ -1,550 +1,686 @@ /*=================================================================== 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 +#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; } +mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromShImage(mitk::ShImage::Pointer mitkImage) +{ + mitk::OdfImage::ItkOdfImageType::Pointer output; + switch (mitkImage->ShOrder()) + { + case 2: + { + typedef itk::ShToOdfImageFilter< float, 2 > ShConverterType; + typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); + mitk::CastToItkImage(mitkImage, itkvol); + typename ShConverterType::Pointer converter = ShConverterType::New(); + converter->SetInput(itkvol); + converter->Update(); + output = converter->GetOutput(); + break; + } + case 4: + { + typedef itk::ShToOdfImageFilter< float, 4 > ShConverterType; + typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); + mitk::CastToItkImage(mitkImage, itkvol); + typename ShConverterType::Pointer converter = ShConverterType::New(); + converter->SetInput(itkvol); + converter->Update(); + output = converter->GetOutput(); + break; + } + case 6: + { + typedef itk::ShToOdfImageFilter< float, 6 > ShConverterType; + typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); + mitk::CastToItkImage(mitkImage, itkvol); + typename ShConverterType::Pointer converter = ShConverterType::New(); + converter->SetInput(itkvol); + converter->Update(); + output = converter->GetOutput(); + break; + } + case 8: + { + typedef itk::ShToOdfImageFilter< float, 8 > ShConverterType; + typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); + mitk::CastToItkImage(mitkImage, itkvol); + typename ShConverterType::Pointer converter = ShConverterType::New(); + converter->SetInput(itkvol); + converter->Update(); + output = converter->GetOutput(); + break; + } + case 10: + { + typedef itk::ShToOdfImageFilter< float, 10 > ShConverterType; + typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); + mitk::CastToItkImage(mitkImage, itkvol); + typename ShConverterType::Pointer converter = ShConverterType::New(); + converter->SetInput(itkvol); + converter->Update(); + output = converter->GetOutput(); + break; + } + case 12: + { + typedef itk::ShToOdfImageFilter< float, 12 > ShConverterType; + typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); + mitk::CastToItkImage(mitkImage, itkvol); + typename ShConverterType::Pointer converter = ShConverterType::New(); + converter->SetInput(itkvol); + converter->Update(); + output = converter->GetOutput(); + break; + } + default: + mitkThrow() << "SH orders higher than 12 are not supported!"; + } + + return output; +} + +mitk::OdfImage::Pointer mitk::convert::GetOdfFromShImage(mitk::ShImage::Pointer mitkImage) +{ + mitk::OdfImage::Pointer image = mitk::OdfImage::New(); + auto img = GetItkOdfFromShImage(mitkImage); + image->InitializeByItk( img.GetPointer() ); + image->SetVolume( img->GetBufferPointer() ); + return image; +} + +mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromTensorImage(mitk::Image::Pointer mitkImage) +{ + typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput( GetItkTensorFromTensorImage(mitkImage) ); + filter->Update(); + return filter->GetOutput(); +} + +mitk::TensorImage::ItkTensorImageType::Pointer mitk::convert::GetItkTensorFromTensorImage(mitk::Image::Pointer mitkImage) +{ + typedef mitk::ImageToItk< mitk::TensorImage::ItkTensorImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(mitkImage); + caster->Update(); + return caster->GetOutput(); +} + +mitk::PeakImage::ItkPeakImageType::Pointer mitk::convert::GetItkPeakFromPeakImage(mitk::Image::Pointer mitkImage) +{ + typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(mitkImage); + caster->Update(); + return caster->GetOutput(); +} + +mitk::OdfImage::Pointer mitk::convert::GetOdfFromTensorImage(mitk::Image::Pointer mitkImage) +{ + mitk::OdfImage::Pointer image = mitk::OdfImage::New(); + auto img = GetItkOdfFromTensorImage(mitkImage); + image->InitializeByItk( img.GetPointer() ); + image->SetVolume( img->GetBufferPointer() ); + return image; +} + +mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromOdfImage(mitk::OdfImage::Pointer mitkImage) +{ + typedef mitk::ImageToItk< mitk::OdfImage::ItkOdfImageType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(mitkImage); + caster->Update(); + return caster->GetOutput(); +} + 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/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h index 79768ec836..593c02520e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h @@ -1,88 +1,87 @@ /*=================================================================== 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 _TrackingHandlerOdf #define _TrackingHandlerOdf #include "mitkTrackingDataHandler.h" #include #include #include #include namespace mitk { /** -* \brief Enables streamline tracking on tensor images. Supports multi tensor tracking by adding multiple tensor images. */ +* \brief Enables streamline tracking on sampled ODF images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerOdf : public TrackingDataHandler { public: TrackingHandlerOdf(); ~TrackingHandlerOdf() override; - typedef TensorImage::PixelType TensorType; typedef OdfImage::ItkOdfImageType ItkOdfImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; void SetSharpenOdfs(bool doSharpen) { m_SharpenOdfs=doSharpen; } void SetOdfThreshold(float odfThreshold){ m_OdfThreshold = odfThreshold; } void SetGfaThreshold(float gfaThreshold){ m_GfaThreshold = gfaThreshold; } void SetOdfImage( ItkOdfImageType::Pointer img ){ m_OdfImage = img; DataModified(); } void SetGfaImage( ItkFloatImgType::Pointer img ){ m_GfaImage = img; DataModified(); } void SetMode( MODE m ) override{ m_Mode = m; } ItkUcharImgType::SpacingType GetSpacing() override{ return m_OdfImage->GetSpacing(); } itk::Point GetOrigin() override{ return m_OdfImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection() override{ return m_OdfImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion() override{ return m_OdfImage->GetLargestPossibleRegion(); } int OdfPower() const; void SetNumProbSamples(int NumProbSamples); bool GetIsOdfFromTensor() const; void SetIsOdfFromTensor(bool OdfFromTensor); protected: int SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles); float m_GfaThreshold; float m_OdfThreshold; bool m_SharpenOdfs; ItkFloatImgType::Pointer m_GfaImage; ///< GFA image used to determine streamline termination. ItkOdfImageType::Pointer m_OdfImage; ///< Input odf image. ItkOdfImageType::Pointer m_WorkingOdfImage; ///< Modified odf image. std::vector< int > m_OdfHemisphereIndices; vnl_matrix< float > m_OdfFloatDirs; int m_NumProbSamples; bool m_OdfFromTensor; itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::Pointer m_GfaInterpolator; itk::LinearInterpolateImageFunction< itk::Image< ItkOdfImageType::PixelType, 3 >, float >::Pointer m_OdfInterpolator; }; } #endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp index ee60e32168..ad77193b5e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp @@ -1,429 +1,425 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkFiberQuantificationView.h" // Qt #include // MITK #include #include #include #include #include #include #include #include // ITK #include #include #include #include #include const std::string QmitkFiberQuantificationView::VIEW_ID = "org.mitk.views.fiberquantification"; using namespace mitk; QmitkFiberQuantificationView::QmitkFiberQuantificationView() : QmitkAbstractView() , m_Controls( 0 ) , m_UpsamplingFactor(5) , m_Visible(false) { } // Destructor QmitkFiberQuantificationView::~QmitkFiberQuantificationView() { } void QmitkFiberQuantificationView::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::QmitkFiberQuantificationViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) ); connect( m_Controls->m_ExtractFiberPeaks, SIGNAL(clicked()), this, SLOT(CalculateFiberDirections()) ); m_Controls->m_TractBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isFib = mitk::TNodePredicateDataType::New(); m_Controls->m_TractBox->SetPredicate( isFib ); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ImageBox->SetZeroEntryText("--"); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3); m_Controls->m_ImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, is3D) ); connect( (QObject*)(m_Controls->m_TractBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect( (QObject*)(m_Controls->m_ImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); } } void QmitkFiberQuantificationView::Activated() { } void QmitkFiberQuantificationView::Deactivated() { } void QmitkFiberQuantificationView::Visible() { m_Visible = true; } void QmitkFiberQuantificationView::Hidden() { m_Visible = false; } void QmitkFiberQuantificationView::SetFocus() { m_Controls->m_ProcessFiberBundleButton->setFocus(); } void QmitkFiberQuantificationView::CalculateFiberDirections() { typedef itk::Image ItkUcharImgType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(m_SelectedFB.back()->GetData()); itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); if (m_SelectedImage.IsNotNull()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(m_SelectedImage, itkMaskImage); fOdfFilter->SetMaskImage(itkMaskImage); } // extract directions from fiber bundle fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetAngularThreshold(cos(m_Controls->m_AngularThreshold->value()*itk::Math::pi/180)); switch (m_Controls->m_FiberDirNormBox->currentIndex()) { case 0: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::GLOBAL_MAX); break; case 1: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::SINGLE_VEC_NORM); break; case 2: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::MAX_VEC_NORM); break; } fOdfFilter->SetSizeThreshold(m_Controls->m_PeakThreshold->value()); fOdfFilter->SetMaxNumDirections(m_Controls->m_MaxNumDirections->value()); fOdfFilter->Update(); QString name = m_SelectedFB.back()->GetName().c_str(); if (m_Controls->m_NumDirectionsBox->isChecked()) { mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( fOdfFilter->GetNumDirectionsImage().GetPointer() ); mitkImage->SetVolume( fOdfFilter->GetNumDirectionsImage()->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(mitkImage); node->SetName((name+"_NUM_DIRECTIONS").toStdString().c_str()); GetDataStorage()->Add(node, m_SelectedFB.back()); } Image::Pointer mitkImage = dynamic_cast(PeakImage::New().GetPointer()); mitk::CastToMitkImage(fOdfFilter->GetDirectionImage(), mitkImage); mitkImage->SetVolume(fOdfFilter->GetDirectionImage()->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(mitkImage); node->SetName( (name+"_DIRECTIONS").toStdString().c_str()); GetDataStorage()->Add(node, m_SelectedFB.back()); } void QmitkFiberQuantificationView::UpdateGui() { m_SelectedFB.clear(); if (m_Controls->m_TractBox->GetSelectedNode().IsNotNull()) m_SelectedFB.push_back(m_Controls->m_TractBox->GetSelectedNode()); m_SelectedImage = nullptr; if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) m_SelectedImage = dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); m_Controls->m_ProcessFiberBundleButton->setEnabled(!m_SelectedFB.empty()); m_Controls->m_ExtractFiberPeaks->setEnabled(!m_SelectedFB.empty()); } void QmitkFiberQuantificationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { UpdateGui(); } void QmitkFiberQuantificationView::ProcessSelectedBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); MITK_WARN("QmitkFiberQuantificationView") << "no fibe bundle selected"; return; } int generationMethod = m_Controls->m_GenerationBox->currentIndex(); for( unsigned int i=0; i(node->GetData())) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); DataNode::Pointer newNode = nullptr; switch(generationMethod){ case 0: newNode = GenerateTractDensityImage(fib, false, true); name += "_TDI"; break; case 1: newNode = GenerateTractDensityImage(fib, false, false); name += "_TDI"; break; case 2: newNode = GenerateTractDensityImage(fib, true, false); name += "_envelope"; break; case 3: newNode = GenerateColorHeatmap(fib); break; case 4: newNode = GenerateFiberEndingsImage(fib); name += "_fiber_endings"; break; case 5: newNode = GenerateFiberEndingsPointSet(fib); name += "_fiber_endings"; break; } if (newNode.IsNotNull()) { newNode->SetName(name.toStdString()); GetDataStorage()->Add(newNode); } } } } // generate pointset displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); int count = 0; int numFibers = fib->GetNumFibers(); for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints>0) { double* point = points->GetPoint(0); itk::Point itkPoint = mitk::imv::GetItkPoint(point); pointSet->InsertPoint(count, itkPoint); count++; } if (numPoints>2) { double* point = points->GetPoint(numPoints-1); itk::Point itkPoint = mitk::imv::GetItkPoint(point); pointSet->InsertPoint(count, itkPoint); count++; } } mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( pointSet ); return node; } // generate image displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib) { typedef unsigned int OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image OutImageType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate rgba heatmap from fiber bundle mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateColorHeatmap(mitk::FiberBundle::Pointer fib) { typedef itk::RGBAPixel OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { itk::Image::Pointer itkImage = itk::Image::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate tract density image from fiber bundle mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute) { mitk::DataNode::Pointer node = mitk::DataNode::New(); if (binary) { typedef unsigned char OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_SelectedImage.IsNotNull()) { mitk::LabelSetImage::Pointer multilabelImage = mitk::LabelSetImage::New(); + multilabelImage->Initialize(m_SelectedImage); multilabelImage->InitializeByLabeledImage(img); - + multilabelImage->GetActiveLabelSet()->SetActiveLabel(1); mitk::Label::Pointer label = multilabelImage->GetActiveLabel(); label->SetName("Tractogram"); -// label->SetColor(color); - label->SetValue(1); -// multilabelImage->GetActiveLabelSet()->AddLabel(label); - multilabelImage->GetActiveLabelSet()->SetActiveLabel(1); - multilabelImage->GetPropertyList()->ConcatenatePropertyList(m_SelectedImage->GetPropertyList()); // init data node node->SetData(multilabelImage); } else { // init data node node->SetData(img); } } else { typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } //generator->SetDoFiberResampling(false); generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node node->SetData(img); } return node; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp index 4f4c4bfce4..d314cba6c4 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.cpp @@ -1,847 +1,808 @@ /*=================================================================== 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. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "QmitkQBallReconstructionView.h" // qt includes #include // itk includes #include "itkTimeProbe.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "itkDiffusionQballReconstructionImageFilter.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include "itkDiffusionMultiShellQballReconstructionImageFilter.h" #include "itkVectorContainer.h" #include "itkB0ImageExtractionImageFilter.h" #include #include "mitkOdfImage.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include #include "mitkDiffusionImagingConfigure.h" #include #include "berryIStructuredSelection.h" #include "berryIWorkbenchWindow.h" #include "berryISelectionService.h" #include #include #include #include +#include const std::string QmitkQBallReconstructionView::VIEW_ID = "org.mitk.views.qballreconstruction"; typedef float TTensorPixelType; const int QmitkQBallReconstructionView::nrconvkernels = 252; struct QbrShellSelection { typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; typedef mitk::DiffusionPropertyHelper::BValueMapType BValueMapType; typedef itk::VectorImage< DiffusionPixelType, 3 > ITKDiffusionImageType; QmitkQBallReconstructionView* m_View; mitk::DataNode * m_Node; std::string m_NodeName; std::vector m_CheckBoxes; QLabel * m_Label; mitk::Image * m_Image; QbrShellSelection(QmitkQBallReconstructionView* view, mitk::DataNode * node) : m_View(view), m_Node(node), m_NodeName(node->GetName()) { m_Image = dynamic_cast (node->GetData()); if(!m_Image) { MITK_ERROR << "QmitkQBallReconstructionView::QbrShellSelection : no image selected"; return; } bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_Node->GetData())) ); if( !isDiffusionImage ) { MITK_ERROR << "QmitkQBallReconstructionView::QbrShellSelection : selected image contains no diffusion information"; return; } GenerateCheckboxes(); } void GenerateCheckboxes() { BValueMapType origMap = mitk::DiffusionPropertyHelper::GetBValueMap(m_Image); BValueMapType::iterator itStart = origMap.begin(); itStart++; BValueMapType::iterator itEnd = origMap.end(); m_Label = new QLabel(m_NodeName.c_str()); m_Label->setVisible(true); m_View->m_Controls->m_QBallSelectionBox->layout()->addWidget(m_Label); for(BValueMapType::iterator it = itStart ; it!= itEnd; it++) { QCheckBox * box = new QCheckBox(QString::number(it->first)); m_View->m_Controls->m_QBallSelectionBox->layout()->addWidget(box); box->setChecked(true); box->setCheckable(true); // box->setVisible(true); m_CheckBoxes.push_back(box); } } void SetVisible(bool vis) { foreach(QCheckBox * box, m_CheckBoxes) { box->setVisible(vis); } } BValueMapType GetBValueSelctionMap() { BValueMapType inputMap = mitk::DiffusionPropertyHelper::GetBValueMap(m_Image); BValueMapType outputMap; unsigned int val = 0; if(inputMap.find(0) == inputMap.end()){ return outputMap; }else{ outputMap[val] = inputMap[val]; } foreach(QCheckBox * box, m_CheckBoxes) { if(box->isChecked()){ val = box->text().toUInt(); outputMap[val] = inputMap[val]; MITK_INFO << val; } } return outputMap; } ~QbrShellSelection() { m_View->m_Controls->m_QBallSelectionBox->layout()->removeWidget(m_Label); delete m_Label; for(std::vector::iterator it = m_CheckBoxes.begin() ; it!= m_CheckBoxes.end(); it++) { m_View->m_Controls->m_QBallSelectionBox->layout()->removeWidget((*it)); delete (*it); } m_CheckBoxes.clear(); } }; QmitkQBallReconstructionView::QmitkQBallReconstructionView() : QmitkAbstractView(), m_Controls(nullptr) { } QmitkQBallReconstructionView::~QmitkQBallReconstructionView() { } void QmitkQBallReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkQBallReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); QStringList items; items << "2" << "4" << "6" << "8" << "10" << "12"; m_Controls->m_QBallReconstructionMaxLLevelComboBox->addItems(items); m_Controls->m_QBallReconstructionMaxLLevelComboBox->setCurrentIndex(1); MethodChoosen(m_Controls->m_QBallReconstructionMethodComboBox->currentIndex()); } } void QmitkQBallReconstructionView::SetFocus() { } void QmitkQBallReconstructionView::CreateConnections() { if ( m_Controls ) { connect( static_cast(m_Controls->m_ButtonStandard), SIGNAL(clicked()), this, SLOT(ReconstructStandard()) ); connect( static_cast(m_Controls->m_QBallReconstructionMethodComboBox), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) ); connect( static_cast(m_Controls->m_QBallReconstructionThreasholdEdit), SIGNAL(valueChanged(int)), this, SLOT(PreviewThreshold(int)) ); connect( static_cast(m_Controls->m_ConvertButton), SIGNAL(clicked()), this, SLOT(ConvertShImage()) ); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New(); m_Controls->m_ImageBox->SetPredicate( isDwi ); m_Controls->m_ShImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isSh = mitk::TNodePredicateDataType::New(); m_Controls->m_ShImageBox->SetPredicate( isSh ); connect( static_cast(m_Controls->m_ImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect( static_cast(m_Controls->m_ShImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); UpdateGui(); } } -template -void QmitkQBallReconstructionView::TemplatedConvertShImage(mitk::ShImage::Pointer mitkImage) -{ - typedef itk::ShToOdfImageFilter< float, ShOrder > ShConverterType; - - typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); - mitk::CastToItkImage(mitkImage, itkvol); - - typename ShConverterType::Pointer converter = ShConverterType::New(); - converter->SetInput(itkvol); - converter->Update(); - - mitk::OdfImage::Pointer image = mitk::OdfImage::New(); - image->InitializeByItk( converter->GetOutput() ); - image->SetVolume( converter->GetOutput()->GetBufferPointer() ); - mitk::DataNode::Pointer node=mitk::DataNode::New(); - node->SetData( image ); - node->SetName(m_Controls->m_ShImageBox->GetSelectedNode()->GetName()); - - GetDataStorage()->Add(node, m_Controls->m_ShImageBox->GetSelectedNode()); -} - void QmitkQBallReconstructionView::ConvertShImage() { if (m_Controls->m_ShImageBox->GetSelectedNode().IsNotNull()) { mitk::ShImage::Pointer mitkImg = dynamic_cast(m_Controls->m_ShImageBox->GetSelectedNode()->GetData()); - switch (mitkImg->ShOrder()) - { - case 2: - TemplatedConvertShImage<2>(mitkImg); - break; - case 4: - TemplatedConvertShImage<4>(mitkImg); - break; - case 6: - TemplatedConvertShImage<6>(mitkImg); - break; - case 8: - TemplatedConvertShImage<8>(mitkImg); - break; - case 10: - TemplatedConvertShImage<10>(mitkImg); - break; - case 12: - TemplatedConvertShImage<12>(mitkImg); - break; - default : - QMessageBox::warning(nullptr, "Error", "Only spherical harmonics orders 2-12 are supported.", QMessageBox::Ok); - } + auto img = mitk::convert::GetOdfFromShImage(mitkImg); + mitk::DataNode::Pointer node= mitk::DataNode::New(); + node->SetData( img ); + node->SetName(m_Controls->m_ShImageBox->GetSelectedNode()->GetName()); + GetDataStorage()->Add(node, m_Controls->m_ShImageBox->GetSelectedNode()); } } void QmitkQBallReconstructionView::UpdateGui() { m_Controls->m_ButtonStandard->setEnabled(false); if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) { m_Controls->m_ButtonStandard->setEnabled(true); GenerateShellSelectionUI(m_Controls->m_ImageBox->GetSelectedNode()); } m_Controls->m_ConvertButton->setEnabled(m_Controls->m_ShImageBox->GetSelectedNode().IsNotNull()); } void QmitkQBallReconstructionView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& /*nodes*/) { UpdateGui(); } void QmitkQBallReconstructionView::Activated() { } void QmitkQBallReconstructionView::Deactivated() { } void QmitkQBallReconstructionView::Visible() { } void QmitkQBallReconstructionView::Hidden() { } void QmitkQBallReconstructionView::ReconstructStandard() { int index = m_Controls->m_QBallReconstructionMethodComboBox->currentIndex(); switch(index) { case 0: { // Numerical Reconstruct(0,0); break; } case 1: { // Standard Reconstruct(1,0); break; } case 2: { // Solid Angle Reconstruct(1,6); break; } case 3: { // ADC Reconstruct(1,4); break; } case 4: { // Raw Signal Reconstruct(1,5); break; } case 5: { // Q-Ball reconstruction Reconstruct(2,0); break; } } } void QmitkQBallReconstructionView::MethodChoosen(int method) { m_Controls->m_QBallSelectionBox->setHidden(true); m_Controls->m_OutputCoeffsImage->setHidden(true); if (method==0) m_Controls->m_ShFrame->setVisible(false); else m_Controls->m_ShFrame->setVisible(true); switch(method) { case 0: m_Controls->m_Description->setText("Numerical recon. (Tuch 2004)"); break; case 1: m_Controls->m_Description->setText("Spherical harmonics recon. (Descoteaux 2007)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 2: m_Controls->m_Description->setText("SH recon. with solid angle consideration (Aganj 2009)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 3: m_Controls->m_Description->setText("SH solid angle with non-neg. constraint (Goh 2009)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 4: m_Controls->m_Description->setText("SH recon. of the plain ADC-profiles"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 5: m_Controls->m_Description->setText("SH recon. of the raw diffusion signal"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 6: m_Controls->m_Description->setText("SH recon. of the multi shell diffusion signal (Aganj 2010)"); m_Controls->m_QBallSelectionBox->setHidden(false); m_Controls->m_OutputCoeffsImage->setHidden(false); break; } } void QmitkQBallReconstructionView::Reconstruct(int method, int normalization) { if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) { if(method == 0) { NumericalQBallReconstruction(m_Controls->m_ImageBox->GetSelectedNode(), normalization); } else if(method == 1) { AnalyticalQBallReconstruction(m_Controls->m_ImageBox->GetSelectedNode(), normalization); } else if(method == 2) { MultiQBallReconstruction(m_Controls->m_ImageBox->GetSelectedNode()); } } } void QmitkQBallReconstructionView::NumericalQBallReconstruction(mitk::DataNode::Pointer node, int normalization) { try { mitk::Image* vols = static_cast(node->GetData()); std::string nodename = node->GetName(); typedef itk::DiffusionQballReconstructionImageFilter QballReconstructionImageFilterType; ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); QballReconstructionImageFilterType::Pointer filter = QballReconstructionImageFilterType::New(); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(vols)); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(vols), itkVectorImagePointer); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); std::string nodePostfix; switch(normalization) { case 0: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); nodePostfix = "_Numerical_Qball"; break; } case 1: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO_B_VALUE); nodePostfix = "_Numerical_ZeroBvalueNormalization_Qball"; break; } case 2: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO); nodePostfix = "_NumericalQball_ZeroNormalization_Qball"; break; } case 3: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_NONE); nodePostfix = "_NumericalQball_NoNormalization_Qball"; break; } default: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); nodePostfix = "_NumericalQball_Qball"; } } filter->Update(); // ODFs TO DATATREE mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer new_node = mitk::DataNode::New(); new_node->SetData( image ); new_node->SetName(nodename+nodePostfix); mitk::ProgressBar::GetInstance()->Progress(); GetDataStorage()->Add(new_node, node); this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return ; } } void QmitkQBallReconstructionView::AnalyticalQBallReconstruction( mitk::DataNode::Pointer node, int normalization) { try { float lambda = m_Controls->m_QBallReconstructionLambdaLineEdit->value(); switch(m_Controls->m_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedAnalyticalQBallReconstruction<2>(node, lambda, normalization); break; } case 1: { TemplatedAnalyticalQBallReconstruction<4>(node, lambda, normalization); break; } case 2: { TemplatedAnalyticalQBallReconstruction<6>(node, lambda, normalization); break; } case 3: { TemplatedAnalyticalQBallReconstruction<8>(node, lambda, normalization); break; } case 4: { TemplatedAnalyticalQBallReconstruction<10>(node, lambda, normalization); break; } case 5: { TemplatedAnalyticalQBallReconstruction<12>(node, lambda, normalization); break; } } this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return; } } template void QmitkQBallReconstructionView::TemplatedAnalyticalQBallReconstruction(mitk::DataNode* dataNodePointer, float lambda, int normalization) { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); mitk::Image* vols = dynamic_cast(dataNodePointer->GetData()); ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(vols)); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(vols), itkVectorImagePointer); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); filter->SetLambda(lambda); std::string nodePostfix; switch(normalization) { case 0: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); nodePostfix = "_SH_Qball"; break; } case 1: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO_B_VALUE); nodePostfix = "_SH_1_Qball"; break; } case 2: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO); nodePostfix = "_SH_2_Qball"; break; } case 3: { filter->SetNormalizationMethod(FilterType::QBAR_NONE); nodePostfix = "_SH_3_Qball"; break; } case 4: { filter->SetNormalizationMethod(FilterType::QBAR_ADC_ONLY); nodePostfix = "_AdcProfile"; break; } case 5: { filter->SetNormalizationMethod(FilterType::QBAR_RAW_SIGNAL); nodePostfix = "_RawSignal"; break; } case 6: { filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); nodePostfix = "_SH_CSA_Qball"; break; } case 7: { filter->SetNormalizationMethod(FilterType::QBAR_NONNEG_SOLID_ANGLE); nodePostfix = "_SH_NonNegCSA_Qball"; break; } default: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); } } filter->Update(); // ODFs TO DATATREE mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetName(dataNodePointer->GetName()+nodePostfix); GetDataStorage()->Add(node, dataNodePointer); if(m_Controls->m_OutputCoeffsImage->isChecked()) { mitk::Image::Pointer coeffsImage = dynamic_cast(mitk::ShImage::New().GetPointer()); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); mitk::DataNode::Pointer coeffsNode=mitk::DataNode::New(); coeffsNode->SetData( coeffsImage ); coeffsNode->SetProperty( "name", mitk::StringProperty::New(dataNodePointer->GetName()+"_SH-Coeffs") ); GetDataStorage()->Add(coeffsNode, node); } } void QmitkQBallReconstructionView::MultiQBallReconstruction(mitk::DataNode::Pointer node) { try { float lambda = m_Controls->m_QBallReconstructionLambdaLineEdit->value(); switch(m_Controls->m_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedMultiQBallReconstruction<2>(lambda, node); break; } case 1: { TemplatedMultiQBallReconstruction<4>(lambda, node); break; } case 2: { TemplatedMultiQBallReconstruction<6>(lambda, node); break; } case 3: { TemplatedMultiQBallReconstruction<8>(lambda, node); break; } case 4: { TemplatedMultiQBallReconstruction<10>(lambda, node); break; } case 5: { TemplatedMultiQBallReconstruction<12>(lambda, node); break; } } } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return ; } } template void QmitkQBallReconstructionView::TemplatedMultiQBallReconstruction(float lambda, mitk::DataNode* dataNodePointer) { typedef itk::DiffusionMultiShellQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); std::string nodename; dataNodePointer->GetStringProperty("name",nodename); mitk::Image* dwi = dynamic_cast(dataNodePointer->GetData()); BValueMapType currSelectionMap = m_ShellSelectorMap[dataNodePointer]->GetBValueSelctionMap(); if(currSelectionMap.size() != 4)// || currSelectionMap.find(0) != currSelectionMap.end()) { QMessageBox::information(0, "Reconstruction not possible:" ,QString("Only three equidistant shells are supported. (ImageName: " + QString(nodename.c_str()) + ")")); return; } BValueMapType::reverse_iterator it1 = currSelectionMap.rbegin(); BValueMapType::reverse_iterator it2 = currSelectionMap.rbegin(); ++it2; // Get average distance int avdistance = 0; for(; it2 != currSelectionMap.rend(); ++it1,++it2) avdistance += (int)it1->first - (int)it2->first; avdistance /= currSelectionMap.size()-1; // Check if all shells are using the same averae distance it1 = currSelectionMap.rbegin(); it2 = currSelectionMap.rbegin(); ++it2; for(; it2 != currSelectionMap.rend(); ++it1,++it2) { if(avdistance != (int)it1->first - (int)it2->first) { QMessageBox::information(0, "Reconstruction not possible:" ,QString("Selected Shells are not in a equidistant configuration. (ImageName: " + QString(nodename.c_str()) + ")")); return; } } ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(dwi, itkVectorImagePointer); filter->SetBValueMap(m_ShellSelectorMap[dataNodePointer]->GetBValueSelctionMap()); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer, mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi)); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); filter->SetLambda(lambda); filter->Update(); // ODFs TO DATATREE mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetName(nodename+"_SH_MultiShell_Qball"); GetDataStorage()->Add(node, dataNodePointer); if(m_Controls->m_OutputCoeffsImage->isChecked()) { mitk::Image::Pointer coeffsImage = dynamic_cast(mitk::ShImage::New().GetPointer()); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); mitk::DataNode::Pointer coeffsNode=mitk::DataNode::New(); coeffsNode->SetData( coeffsImage ); coeffsNode->SetProperty( "name", mitk::StringProperty::New(dataNodePointer->GetName()+"_SH-Coeffs") ); GetDataStorage()->Add(coeffsNode, node); } } void QmitkQBallReconstructionView::GenerateShellSelectionUI(mitk::DataNode::Pointer node) { std::map tempMap; if(m_ShellSelectorMap.find( node.GetPointer() ) != m_ShellSelectorMap.end()) { tempMap[node.GetPointer()] = m_ShellSelectorMap[node.GetPointer()]; m_ShellSelectorMap.erase(node.GetPointer()); } else { tempMap[node.GetPointer()] = new QbrShellSelection(this, node ); tempMap[node.GetPointer()]->SetVisible(true); } for(std::map::iterator it = m_ShellSelectorMap.begin(); it != m_ShellSelectorMap.end();it ++) { delete it->second; } m_ShellSelectorMap.clear(); m_ShellSelectorMap = tempMap; } void QmitkQBallReconstructionView::PreviewThreshold(int threshold) { if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) { mitk::Image* vols = static_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); // Extract b0 image ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filterB0 = FilterType::New(); filterB0->SetInput( itkVectorImagePointer ); filterB0->SetDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(vols)); filterB0->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); typedef itk::Image ImageType; typedef itk::Image SegmentationType; typedef itk::BinaryThresholdImageFilter ThresholdFilterType; // apply threshold ThresholdFilterType::Pointer filterThreshold = ThresholdFilterType::New(); filterThreshold->SetInput(filterB0->GetOutput()); filterThreshold->SetLowerThreshold(threshold); filterThreshold->SetInsideValue(0); filterThreshold->SetOutsideValue(1); // mark cut off values red filterThreshold->Update(); mitkImage->InitializeByItk( filterThreshold->GetOutput() ); mitkImage->SetVolume( filterThreshold->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node; if (this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", m_Controls->m_ImageBox->GetSelectedNode())) { node = this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", m_Controls->m_ImageBox->GetSelectedNode()); } else { // create a new node, to show thresholded values node = mitk::DataNode::New(); GetDataStorage()->Add( node, m_Controls->m_ImageBox->GetSelectedNode() ); node->SetProperty( "name", mitk::StringProperty::New("ThresholdOverlay")); node->SetBoolProperty("helper object", true); } node->SetData( mitkImage ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.h index e8a73bebc4..8cf890cd0b 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkQBallReconstructionView.h @@ -1,136 +1,133 @@ /*=================================================================== 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 _QMITKQBALLRECONSTRUCTIONVIEW_H_INCLUDED #define _QMITKQBALLRECONSTRUCTIONVIEW_H_INCLUDED #include #include #include #include "ui_QmitkQBallReconstructionViewControls.h" #include #include #include #include typedef short DiffusionPixelType; struct QbrSelListener; struct QbrShellSelection; /*! * \ingroup org_mitk_gui_qt_qballreconstruction_internal * * \brief QmitkQBallReconstructionView * * Document your class here. */ class QmitkQBallReconstructionView : public QmitkAbstractView, public mitk::ILifecycleAwarePart { friend struct QbrSelListener; friend struct QbrShellSelection; typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; typedef mitk::DiffusionPropertyHelper::BValueMapType BValueMapType; typedef itk::VectorImage< DiffusionPixelType, 3 > ITKDiffusionImageType; // this is needed for all Qt objects that should have a MOC object (everything that derives from QObject) Q_OBJECT public: static const std::string VIEW_ID; QmitkQBallReconstructionView(); virtual ~QmitkQBallReconstructionView(); virtual void CreateQtPartControl(QWidget *parent) override; /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; /// \brief Called when the view gets activated virtual void Activated() override; /// \brief Called when the view gets deactivated virtual void Deactivated() override; /// \brief Called when the view becomes visible virtual void Visible() override; /// \brief Called when the view becomes hidden virtual void Hidden() override; static const int nrconvkernels; protected slots: void UpdateGui(); void ReconstructStandard(); void ConvertShImage(); void MethodChoosen(int method); void Reconstruct(int method, int normalization); void NumericalQBallReconstruction(mitk::DataNode::Pointer node, int normalization); void AnalyticalQBallReconstruction(mitk::DataNode::Pointer node, int normalization); void MultiQBallReconstruction(mitk::DataNode::Pointer node); /** * @brief PreviewThreshold Generates a preview of the values that are cut off by the thresholds * @param threshold */ void PreviewThreshold(int threshold); protected: - template - void TemplatedConvertShImage(mitk::ShImage::Pointer mitkImage); - /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkQBallReconstructionViewControls* m_Controls; template void TemplatedAnalyticalQBallReconstruction(mitk::DataNode* dataNodePointer, float lambda, int normalization); template void TemplatedMultiQBallReconstruction(float lambda, mitk::DataNode*); private: std::map< const mitk::DataNode *, QbrShellSelection * > m_ShellSelectorMap; void GenerateShellSelectionUI(mitk::DataNode::Pointer node); }; #endif // _QMITKQBALLRECONSTRUCTIONVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index 8a25387acb..bfe4a2eb0c 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,1045 +1,1042 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include // VTK #include #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) : m_View(view) { } void QmitkStreamlineTrackingWorker::run() { m_View->m_Tracker->Update(); m_View->m_TrackingThread.quit(); } QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_TrackingWorker(this) , m_Controls(nullptr) , m_FirstTensorProbRun(true) , m_FirstInteractiveRun(true) , m_TrackingHandler(nullptr) , m_ThreadIsRunning(false) , m_DeleteTrackingHandler(false) , m_Visible(false) , m_LastPrior("") , m_TrackingPriorHandler(nullptr) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_TargetImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_PriorImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_ForestSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_ExclusionImageSelectionWidget->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isPeakImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestSelectionWidget->SetNodePredicate(isTractographyForest); m_Controls->m_FaImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_FaImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_SeedImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_SeedImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_SeedImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_MaskImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_MaskImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_MaskImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_StopImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_StopImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_StopImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_TargetImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_TargetImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_TargetImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_PriorImageSelectionWidget->SetNodePredicate( isPeakImagePredicate ); m_Controls->m_PriorImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_PriorImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_ExclusionImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_ExclusionImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_ExclusionImageSelectionWidget->SetSelectionIsOptional(true); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_SeedImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TargetImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_PriorImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ExclusionImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskInterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_LoopCheckBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TrialsPerSeedBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_EpConstraintsBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); m_Controls->m_SeedsPerVoxelBox->editingFinished(); m_Controls->m_NumFibersBox->editingFinished(); m_Controls->m_ScalarThresholdBox->editingFinished(); m_Controls->m_OdfCutoffBox->editingFinished(); m_Controls->m_StepSizeBox->editingFinished(); m_Controls->m_SamplingDistanceBox->editingFinished(); m_Controls->m_AngularThresholdBox->editingFinished(); m_Controls->m_MinTractLengthBox->editingFinished(); m_Controls->m_fBox->editingFinished(); m_Controls->m_gBox->editingFinished(); m_Controls->m_NumSamplesBox->editingFinished(); m_Controls->m_SeedRadiusBox->editingFinished(); m_Controls->m_NumSeedsBox->editingFinished(); m_Controls->m_LoopCheckBox->editingFinished(); m_Controls->m_TrialsPerSeedBox->editingFinished(); StartStopTrackingGui(false); } UpdateGui(); } void QmitkStreamlineTrackingView::StopTractography() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); } void QmitkStreamlineTrackingView::TimerUpdate() { if (m_Tracker.IsNull()) return; QString status_text(m_Tracker->GetStatusText().c_str()); m_Controls->m_StatusTextBox->setText(status_text); } void QmitkStreamlineTrackingView::BeforeThread() { m_TrackingTimer->start(1000); } void QmitkStreamlineTrackingView::AfterThread() { m_TrackingTimer->stop(); if (!m_Tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_ParentNode->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked() && fiberBundle->GetNumberOfLines()>0) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); m_Tracker->SetDicomProperties(fib); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_ParentNode->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); node->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); GetDataStorage()->Add(node, m_ParentNode); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_ParentNode->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_ParentNode); } } if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if(!CheckAndStoreLastParams(sender()) && !posChanged) return; if (m_ThreadIsRunning || !m_Visible) return; if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked()) ) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); DoFiberTracking(); } bool QmitkStreamlineTrackingView::CheckAndStoreLastParams(QObject* obj) { if (obj!=nullptr) { std::string new_val = ""; if(qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else if (qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else return true; if (m_LastTractoParams.find(obj->objectName())==m_LastTractoParams.end()) { m_LastTractoParams[obj->objectName()] = new_val; return false; } else if (m_LastTractoParams.at(obj->objectName()) != new_val) { m_LastTractoParams[obj->objectName()] = new_val; return true; } else if (m_LastTractoParams.at(obj->objectName()) == new_val) return false; } return true; } void QmitkStreamlineTrackingView::OnParameterChanged() { UpdateGui(); if(!CheckAndStoreLastParams(sender())) return; if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageSelectionWidget->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { if (m_FirstInteractiveRun) { QMessageBox::information(nullptr, "Information", "Place and move a spherical seed region anywhere in the image by left-clicking and dragging. If the seed region is colored red, tracking is in progress. If the seed region is colored white, tracking is finished.\nPlacing the seed region for the first time in a newly selected dataset might cause a short delay, since the tracker needs to be initialized."); m_FirstInteractiveRun = false; } QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); disconnect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } } void QmitkStreamlineTrackingView::Activated() { } void QmitkStreamlineTrackingView::Deactivated() { } void QmitkStreamlineTrackingView::Visible() { m_Visible = true; } void QmitkStreamlineTrackingView::Hidden() { m_Visible = false; m_Controls->m_InteractiveBox->setChecked(false); ToggleInteractive(); } void QmitkStreamlineTrackingView::OnSliceChanged() { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (!m_ThreadIsRunning && m_TrackingHandler != nullptr) { if (m_TrackingPriorHandler != nullptr) { delete m_TrackingPriorHandler; m_TrackingPriorHandler = nullptr; } delete m_TrackingHandler; m_TrackingHandler = nullptr; m_DeleteTrackingHandler = false; m_LastPrior = ""; } else if (m_ThreadIsRunning) { m_DeleteTrackingHandler = true; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { std::vector< mitk::DataNode::Pointer > last_nodes = m_InputImageNodes; m_InputImageNodes.clear(); m_InputImages.clear(); m_AdditionalInputImages.clear(); bool retrack = false; for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { - if( dynamic_cast(node->GetData()) ) - { - m_InputImageNodes.push_back(node); - m_InputImages.push_back(dynamic_cast(node->GetData())); - retrack = true; - } - else if ( dynamic_cast(node->GetData()) ) - { - m_InputImageNodes.push_back(node); - m_InputImages.push_back(dynamic_cast(node->GetData())); - retrack = true; - } - else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData())) ) + if( dynamic_cast(node->GetData()) || + dynamic_cast(node->GetData()) || + dynamic_cast(node->GetData()) || + mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData()))) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr) { - int dim = img->GetDimension(); + auto dim = img->GetDimension(); unsigned int* dimensions = img->GetDimensions(); if (dim==4 && dimensions[3]%3==0) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if (dim==3) { m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } } } // sometimes the OnSelectionChanged event is sent twice and actually no selection has changed for the first event. We need to catch that. if (last_nodes.size() == m_InputImageNodes.size()) { bool same_nodes = true; for (unsigned int i=0; im_TensorImageLabel->setText("select in data-manager"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageSelectionWidget->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestSelectionWidget->setVisible(false); m_Controls->m_ForestLabel->setVisible(false); m_Controls->commandLinkButton->setEnabled(false); m_Controls->m_TrialsPerSeedBox->setEnabled(false); m_Controls->m_TrialsPerSeedLabel->setEnabled(false); m_Controls->m_TargetImageSelectionWidget->setEnabled(false); m_Controls->m_TargetImageLabel->setEnabled(false); if (m_Controls->m_InteractiveBox->isChecked()) { m_Controls->m_InteractiveSeedingFrame->setVisible(true); m_Controls->m_StaticSeedingFrame->setVisible(false); m_Controls->commandLinkButton_2->setVisible(false); m_Controls->commandLinkButton->setVisible(false); } else { m_Controls->m_InteractiveSeedingFrame->setVisible(false); m_Controls->m_StaticSeedingFrame->setVisible(true); m_Controls->commandLinkButton_2->setVisible(m_ThreadIsRunning); m_Controls->commandLinkButton->setVisible(!m_ThreadIsRunning); } if (m_Controls->m_EpConstraintsBox->currentIndex()>0) { m_Controls->m_TargetImageSelectionWidget->setEnabled(true); m_Controls->m_TargetImageLabel->setEnabled(true); } // trials per seed are only important for probabilistic tractography if (m_Controls->m_ModeBox->currentIndex()==1) { m_Controls->m_TrialsPerSeedBox->setEnabled(true); m_Controls->m_TrialsPerSeedLabel->setEnabled(true); } if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText( ( std::to_string(m_InputImageNodes.size()) + " images selected").c_str() ); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestSelectionWidget->setVisible(true); m_Controls->m_ForestLabel->setVisible(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } } void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) { m_ThreadIsRunning = start; if (!m_Controls->m_InteractiveBox->isChecked()) { m_Controls->commandLinkButton_2->setVisible(start); m_Controls->commandLinkButton->setVisible(!start); m_Controls->m_InteractiveBox->setEnabled(!start); m_Controls->m_StatusTextBox->setVisible(start); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_InputImages.empty()) { QMessageBox::information(nullptr, "Information", "Please select an input image in the datamaneger (tensor, ODF, peak or dMRI image)!"); return; } if (m_ThreadIsRunning || !m_Visible) return; if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) return; StartStopTrackingGui(true); m_Tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImages.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::TensorImage::ItkTensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkImg ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); dynamic_cast(m_TrackingHandler)->SetIsOdfFromTensor(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (int i=0; i<(int)m_InputImages.size(); i++) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(i)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(m_TrackingHandler)->AddTensorImage(itkImg); } if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(m_TrackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } - else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) + else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) || + dynamic_cast(m_InputImageNodes.at(0)->GetData())) { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); - mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); - mitk::CastToItkImage(m_InputImages.at(0), itkImg); - dynamic_cast(m_TrackingHandler)->SetOdfImage(itkImg); + + if (dynamic_cast(m_InputImageNodes.at(0)->GetData())) + dynamic_cast(m_TrackingHandler)->SetOdfImage(mitk::convert::GetItkOdfFromShImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); + else + { + mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); + mitk::CastToItkImage(m_InputImages.at(0), itkImg); + } if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestSelectionWidget->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography (raw dMRI tractography) selected. Did you accidentally select the raw diffusion-weighted image in the datamanager?"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestSelectionWidget->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { int num_previous_directions = (forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { int num_previous_directions = (forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); StartStopTrackingGui(false); return; } } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); StartStopTrackingGui(false); return; } try { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(0)); caster->SetCopyMemFlag(true); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(itkImg); } dynamic_cast(m_TrackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); } catch(...) { QMessageBox::information(nullptr, "Error", "Peak tracker could not be initialized. Is your input image in the correct format (4D float image, peaks in the 4th dimension)?"); StartStopTrackingGui(false); return; } } m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { m_Tracker->SetSeedPoints(m_SeedPoints); } else if (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetSeedImage(mask); } if (m_Controls->m_MaskImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetMaskImage(mask); } if (m_Controls->m_StopImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetStoppingRegions(mask); } if (m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetTargetRegions(mask); } if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNotNull()) { if (m_LastPrior!=m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetUID() || m_TrackingPriorHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(dynamic_cast(m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetData())); caster->SetCopyMemFlag(true); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingPriorHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingPriorHandler)->SetPeakImage(itkImg); dynamic_cast(m_TrackingPriorHandler)->SetPeakThreshold(0.0); m_LastPrior = m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetUID(); } m_TrackingPriorHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); m_TrackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); m_Tracker->SetTrackingPriorHandler(m_TrackingPriorHandler); m_Tracker->SetTrackingPriorWeight(m_Controls->m_PriorWeightBox->value()); m_Tracker->SetTrackingPriorAsMask(m_Controls->m_PriorAsMaskBox->isChecked()); m_Tracker->SetIntroduceDirectionsFromPrior(m_Controls->m_NewDirectionsFromPriorBox->isChecked()); } else if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNull()) m_Tracker->SetTrackingPriorHandler(nullptr); if (m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetExclusionRegions(mask); } // Endpoint constraints switch (m_Controls->m_EpConstraintsBox->currentIndex()) { case 0: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); m_Tracker->SetTargetRegions(nullptr); break; case 1: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); break; case 2: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); break; case 3: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); break; case 4: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); break; case 5: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); break; case 6: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); break; } if (m_Tracker->GetEndpointConstraint()!=itk::StreamlineTrackingFilter::EndpointConstraints::NONE && m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Error", "Endpoint constraints are used but no target image is set!"); StartStopTrackingGui(false); return; } else if (m_Tracker->GetEndpointConstraint()==itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET && (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNull()|| m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) ) { QMessageBox::information(nullptr, "Error", "Endpoint constraint EPS_IN_SEED_AND_TARGET is used but no target or no seed image is set!"); StartStopTrackingGui(false); return; } m_Tracker->SetInterpolateMasks(m_Controls->m_MaskInterpolationBox->isChecked()); m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); m_Tracker->SetTrialsPerSeed(m_Controls->m_TrialsPerSeedBox->value()); m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); m_Tracker->SetLoopCheck(m_Controls->m_LoopCheckBox->value()); m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); m_ParentNode = m_InputImageNodes.at(0); m_TrackingThread.start(QThread::LowestPriority); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingDataImportPage.dox b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingDataImportPage.dox index 304151e059..d2ef1036fb 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingDataImportPage.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingDataImportPage.dox @@ -1,82 +1,82 @@ /** \page QmitkDiffusionImagingDataImportPage Data formats, import and export MITK Diffusion supports many standard image formats such as NIFTI, NRRD and DICOM as well as common tractography file formats such as vtk/fib, trk, tck, and tractography DICOM. By including multiple tractography file formats native to other tools (e.g. 3D Slicer, MRtrix, DIPY), MITK Diffusion integrates seamlessly in complex workflows involving other tools. Additional file types such as spherical-harmonic coefficient files and voxel-wise fiber orientation or peak images are compliant with MRtrix. Data can be loaded using the open file dialog via the menu bar or by simply dragging and dropping the file into the data manager or one of data display windows. MITK Diffusion uses a left-posterior-superior coordinate system convention. Other toolkits, e.g. MRtrix, use an RAS coordinate system. This can cause flips for example of the complete tractogram or of the voxel-wise peaks or ODFs. MITK Diffusion implements mechanisms to deal with this by automatically catching and converting some cases and by providing the tools to manually correct for it, but it is important to keep this issue in mind and to perform the data processing in a corresponding thorough manner. \section DWI_format Diffusion-weighted Images This section describes the different file formats of raw diffusion-weighted images that are supported by MITK Diffusion. General points: \li MITK Diffusion internally applies the image rotation matrix to the diffusion-gradient vectors while loading the image. The original gradient directions are retained and used when saving the image again. \subsection dwi_nrrd NRRD (.nrrd/.dwi) NRRD (http://teem.sourceforge.net/nrrd/format.html) is the default file format for saving diffusion-weighted images with MITK Diffusion. The gradient and b-value information is directly stored in the file header. The image information is stored as a 3D vector image. Storing the file as .dwi or .nrrd is equivalent, only the file ending is different. Diffusion-weighted images are discerned from other images by the tag "modality:=DWMRI" in the NRRD image header. The gradient information is stored in the NRRD header in the following way: DWMRI_b-value:=1000.000000 DWMRI_gradient_0000:=0.000000 0.000000 0.000000 DWMRI_gradient_0001:=0.000000 0.000000 1.000000 DWMRI_gradient_0002:=-0.051773 0.252917 0.966102 ... The b-value of the individual gradient directions is encoded via the squared norm of the respective vector. The tag "DWMRI_b-value" defines which b-value corresponds to a vector with norm 1. In some cases, dMRI NRRD files contain an additional "measurement frame" matrix, which specifies an additional rotation of the gradient vectors. This matrix is automatically applied when loading the file. The original gradient directions are retained and used when saving the image again. \subsection dwi_nifti NIFTI (.nii/.ni.gz) Diffusion-weighted images can be saved and loaded as NIFTI files (NIFTI-1 https://nifti.nimh.nih.gov/nifti-1). The gradient vector information is stored in to separate files for b-values (filename.bvals) and gradient vectors (filename.bvecs). When loading a nifti file (.nii or .nii.gz), MITK Diffusion looks for these two additional files (.bval/.bvals and .bvec/.bvecs) and if they are found, MITK will offer to load the image as diffusion-weighted image. In contrast to the NRRD format, all gradient vectors should have a length of 1. The b-values are stored explicietly in the .bval/.bvals file. When loaded into MITK, this information is again encoded into the gradient vectors, as it is done in the NRRD file format. \subsection dwi_dicom DICOM MITK is capable of importing diffusion-weighted DICOM images from GE, Siemens and Philips. Writing images is DICOM format is NOT supported! Mosaic images can be directly converted to regular images during import. To load a dMRI DICOM, simply drag and drop any file of the series into the application. \section special_format Special Image Types \subsection DTI_format Diffusion Tensor Images (.dti) The default format for diffusion tensor images in MITK Diffusion is a 3D NRRD file with a "3D-symmetric-matrix" pixel type. A tensor is encoded as 6 float values. The file ending for diffusion tensor files is .dti. MITK Diffusion is also able to read NIFTI DTI files (6 or 9 component format), which are for example generated by the Camino multi tensor reconstruction. To be recognized as DTI files, the .nii or .nii.gz files have to be renamed to .dti. \subsection ODF_format ODF Images (.odf, .qbi (deprecated) ) MITK stores ODFs as 252 float values spherically sampled from the continuous ODF. The sampling directions are generate by a 5-fold subdivisions of an icosahedron. ODF images are stored in NRRD file format with the ending .odf. The image information is stored as a 3D vector image with a vector length of 252. The specific ODF sampling directions can be found here. \subsection sh_format Spherical Harmonics (.nii, .nii.gz, .nrrd) -Many applications in dMRI are using spherical harmonics to store spherical functions such as ODFs. MITK Diffusion stores these files as 4D float images this is the same format as MRtrix is using. SH Images connot be used for tractography directly but they must be converted into sampled ODF images. The conversion tool is located in the \ref org_mitk_views_qballreconstruction view. +Many applications in dMRI are using spherical harmonics to store spherical functions such as ODFs. MITK Diffusion stores these files as 4D float images this is the same format as MRtrix is using. \subsection peak_format Peak Images (.nii, .nii.gz, .nrrd) MITK Diffusion stores peak images, resulting for example from an ODF maxima extraction, as 4D float images. The peak vector components are stored in the 4th dimension, therefore dimension 4 always contains a multiple of 3 entries. This format is the same as the format used by MRtrix. \section Tract_format Tractography Formats MITK Diffusion supports multiple formats for tractography files that are commonly used in the dMRI community. As in all other toolkits, the fiber point coordinates are stored as physical/world coordinates without any additional transformation. \subsection tract_vtk VTK (.vtk/.fib) The default format for tractograms is VTK (vtkPolyData) with the file endings .fib or .vtk. The advantage of the VTK format is that it can store additional information such as fiber weights and fiber colors. By default, both are saved with the actual tract information. \subsection tract_trk TrackVis (.trk) TRK is the tractography file format used by TrackVis and DIPY. See http://www.trackvis.org/docs/?subsect=fileformat for a detailed description of a format. \subsection tract_tck MRtrix (.tck) MITK Diffusion is able to read the tck file format native to MRtrix. Writing this format is currently not supported. By default, MRtrix uses a RAS coordinate convention in contrast to the MITK Diffusion (and also ITK) convention of LPS. To compensate for this, the fiber coordinates are negated in the x and y dimension. \subsection tract_dicom DICOM Tractography DICOM files compliant with the supplement 181 of the DICOM standard can be read and written by MITK Diffusion. This format is for example supported by the neuronavigation software of Brainlab. DICOM tags of the read tractogram can be manually set or modified using the "Properties" view in MITK Diffusion. To do this, enable the "Developer Mode" option in Window>Preferences>Properties>Developer Mode. Then select the fiber bundle in the data manager and select the "Base Data" property list in the corresponding combobox of the "Properties" view. */