diff --git a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/Properties/mitkDiffusionPropertyHelper.cpp b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/Properties/mitkDiffusionPropertyHelper.cpp index c64608bbeb..6d841480f9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/Properties/mitkDiffusionPropertyHelper.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/Properties/mitkDiffusionPropertyHelper.cpp @@ -1,568 +1,568 @@ /*=================================================================== 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 "mitkDiffusionPropertyHelper.h" #include #include #include #include #include const std::string mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME = "DWMRI.GradientDirections"; const std::string mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME = "DWMRI.OriginalGradientDirections"; const std::string mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME = "DWMRI.MeasurementFrame"; const std::string mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME = "DWMRI.ReferenceBValue"; const std::string mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME = "DWMRI.BValueMap"; const std::string mitk::DiffusionPropertyHelper::MODALITY = "DWMRI.Modality"; const std::string mitk::DiffusionPropertyHelper::KEEP_ORIGINAL_DIRECTIONS = "DWMRI.KeepOriginalDirections"; mitk::DiffusionPropertyHelper::DiffusionPropertyHelper() { } mitk::DiffusionPropertyHelper::DiffusionPropertyHelper( mitk::Image* inputImage) : m_Image( inputImage ) { // Update props } mitk::DiffusionPropertyHelper::~DiffusionPropertyHelper() { } mitk::DiffusionPropertyHelper::ImageType::Pointer mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk::Image* image) { ImageType::Pointer vectorImage = ImageType::New(); mitk::CastToItkImage(image, vectorImage); return vectorImage; } mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::DiffusionPropertyHelper::CalcAveragedDirectionSet(double precision, GradientDirectionsContainerType::Pointer directions) { // save old and construct new direction container GradientDirectionsContainerType::Pointer newDirections = GradientDirectionsContainerType::New(); // fill new direction container for(GradientDirectionsContainerType::ConstIterator gdcitOld = directions->Begin(); gdcitOld != directions->End(); ++gdcitOld) { // already exists? bool found = false; for(GradientDirectionsContainerType::ConstIterator gdcitNew = newDirections->Begin(); gdcitNew != newDirections->End(); ++gdcitNew) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { found = true; break; } } // if not found, add it to new container if(!found) { newDirections->push_back(gdcitOld.Value()); } } return newDirections; } void mitk::DiffusionPropertyHelper::AverageRedundantGradients(double precision) { mitk::GradientDirectionsProperty* DirectionsProperty = static_cast( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ); GradientDirectionsContainerType::Pointer oldDirs = DirectionsProperty->GetGradientDirectionsContainer(); GradientDirectionsContainerType::Pointer newDirs = CalcAveragedDirectionSet(precision, oldDirs); // if sizes equal, we do not need to do anything in this function if(oldDirs->size() == newDirs->size()) return; // new image ImageType::Pointer oldImage = ImageType::New(); mitk::CastToItkImage( m_Image, oldImage); ImageType::Pointer newITKImage = ImageType::New(); newITKImage->SetSpacing( oldImage->GetSpacing() ); // Set the image spacing newITKImage->SetOrigin( oldImage->GetOrigin() ); // Set the image origin newITKImage->SetDirection( oldImage->GetDirection() ); // Set the image direction newITKImage->SetLargestPossibleRegion( oldImage->GetLargestPossibleRegion() ); newITKImage->SetVectorLength( newDirs->size() ); newITKImage->SetBufferedRegion( oldImage->GetLargestPossibleRegion() ); newITKImage->Allocate(); // average image data that corresponds to identical directions itk::ImageRegionIterator< ImageType > newIt(newITKImage, newITKImage->GetLargestPossibleRegion()); newIt.GoToBegin(); itk::ImageRegionIterator< ImageType > oldIt(oldImage, oldImage->GetLargestPossibleRegion()); oldIt.GoToBegin(); // initial new value of voxel ImageType::PixelType newVec; newVec.SetSize(newDirs->size()); newVec.AllocateElements(newDirs->size()); // find which gradients should be averaged GradientDirectionsContainerType::Pointer oldDirections = oldDirs; std::vector > dirIndices; for(GradientDirectionsContainerType::ConstIterator gdcitNew = newDirs->Begin(); gdcitNew != newDirs->End(); ++gdcitNew) { dirIndices.push_back(std::vector(0)); for(GradientDirectionsContainerType::ConstIterator gdcitOld = oldDirs->Begin(); gdcitOld != oldDirections->End(); ++gdcitOld) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { //MITK_INFO << gdcitNew.Value() << " " << gdcitOld.Value(); dirIndices[gdcitNew.Index()].push_back(gdcitOld.Index()); } } } //int ind1 = -1; while(!newIt.IsAtEnd()) { // progress //typename ImageType::IndexType ind = newIt.GetIndex(); //ind1 = ind.m_Index[2]; // init new vector with zeros newVec.Fill(0.0); // the old voxel value with duplicates ImageType::PixelType oldVec = oldIt.Get(); for(unsigned int i=0; iGetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( newDirs ) ); m_Image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( newDirs ) ); ApplyMeasurementFrameAndRotationMatrix(); UpdateBValueMap(); std::cout << std::endl; } void mitk::DiffusionPropertyHelper::ApplyMeasurementFrameAndRotationMatrix() { if( m_Image->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).IsNull() ) { return; } GradientDirectionsContainerType::Pointer originalDirections = static_cast( m_Image->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer())->GetGradientDirectionsContainer(); MeasurementFrameType measurementFrame = GetMeasurementFrame(m_Image); mitk::Vector3D s = m_Image->GetGeometry()->GetSpacing(); mitk::VnlVector c0 = m_Image->GetGeometry()->GetMatrixColumn(0)/s[0]; mitk::VnlVector c1 = m_Image->GetGeometry()->GetMatrixColumn(1)/s[1]; mitk::VnlVector c2 = m_Image->GetGeometry()->GetMatrixColumn(2)/s[2]; MeasurementFrameType imageRotationMatrix; imageRotationMatrix[0][0] = c0[0]; imageRotationMatrix[1][0] = c0[1]; imageRotationMatrix[2][0] = c0[2]; imageRotationMatrix[0][1] = c1[0]; imageRotationMatrix[1][1] = c1[1]; imageRotationMatrix[2][1] = c1[2]; imageRotationMatrix[0][2] = c2[0]; imageRotationMatrix[1][2] = c2[1]; imageRotationMatrix[2][2] = c2[2]; GradientDirectionsContainerType::Pointer directions = GradientDirectionsContainerType::New(); if( originalDirections.IsNull() || ( originalDirections->size() == 0 ) ) { // original direction container was not set return; } bool keep_originals = false; m_Image->GetPropertyList()->GetBoolProperty(mitk::DiffusionPropertyHelper::KEEP_ORIGINAL_DIRECTIONS.c_str(), keep_originals); if (!keep_originals) { MITK_INFO << "Applying measurement frame to diffusion-gradient directions:"; std::cout << measurementFrame << std::endl; - MITK_INFO << "Applying image totation to diffusion-gradient directions:"; + MITK_INFO << "Applying image rotation to diffusion-gradient directions:"; std::cout << imageRotationMatrix << std::endl; } int c = 0; for(GradientDirectionsContainerType::ConstIterator gdcit = originalDirections->Begin(); gdcit != originalDirections->End(); ++gdcit) { vnl_vector vec = gdcit.Value(); if (!keep_originals) { vec = vec.pre_multiply(measurementFrame); vec = vec.pre_multiply(imageRotationMatrix); } directions->InsertElement(c, vec); c++; } m_Image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( directions ) ); } void mitk::DiffusionPropertyHelper::UnApplyMeasurementFrameAndRotationMatrix() { if( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).IsNull() ) { return; } GradientDirectionsContainerType::Pointer modifiedDirections = static_cast( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer())->GetGradientDirectionsContainer(); MeasurementFrameType measurementFrame = GetMeasurementFrame(m_Image); measurementFrame = vnl_matrix_inverse(measurementFrame).pinverse(); mitk::Vector3D s = m_Image->GetGeometry()->GetSpacing(); mitk::VnlVector c0 = m_Image->GetGeometry()->GetMatrixColumn(0)/s[0]; mitk::VnlVector c1 = m_Image->GetGeometry()->GetMatrixColumn(1)/s[1]; mitk::VnlVector c2 = m_Image->GetGeometry()->GetMatrixColumn(2)/s[2]; MeasurementFrameType imageRotationMatrix; imageRotationMatrix[0][0] = c0[0]; imageRotationMatrix[1][0] = c0[1]; imageRotationMatrix[2][0] = c0[2]; imageRotationMatrix[0][1] = c1[0]; imageRotationMatrix[1][1] = c1[1]; imageRotationMatrix[2][1] = c1[2]; imageRotationMatrix[0][2] = c2[0]; imageRotationMatrix[1][2] = c2[1]; imageRotationMatrix[2][2] = c2[2]; imageRotationMatrix = vnl_matrix_inverse(imageRotationMatrix).pinverse(); GradientDirectionsContainerType::Pointer directions = GradientDirectionsContainerType::New(); if( modifiedDirections.IsNull() || ( modifiedDirections->size() == 0 ) ) { // original direction container was not set return; } bool keep_originals = false; m_Image->GetPropertyList()->GetBoolProperty(mitk::DiffusionPropertyHelper::KEEP_ORIGINAL_DIRECTIONS.c_str(), keep_originals); if (!keep_originals) { - MITK_INFO << "Reverting image totation to diffusion-gradient directions:"; + MITK_INFO << "Reverting image rotation to diffusion-gradient directions:"; std::cout << imageRotationMatrix << std::endl; MITK_INFO << "Reverting measurement frame to diffusion-gradient directions:"; std::cout << measurementFrame << std::endl; } int c = 0; for(GradientDirectionsContainerType::ConstIterator gdcit = modifiedDirections->Begin(); gdcit != modifiedDirections->End(); ++gdcit) { vnl_vector vec = gdcit.Value(); if (!keep_originals) { vec = vec.pre_multiply(imageRotationMatrix); vec = vec.pre_multiply(measurementFrame); } directions->InsertElement(c, vec); c++; } m_Image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( directions ) ); } void mitk::DiffusionPropertyHelper::ClearMeasurementFrameAndRotationMatrixFromGradients(mitk::Image* image) { if( image->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).IsNull() ) { return; } GradientDirectionsContainerType::Pointer originalDirections = static_cast( image->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer())->GetGradientDirectionsContainer(); GradientDirectionsContainerType::Pointer directions = GradientDirectionsContainerType::New(); if( originalDirections.IsNull() || ( originalDirections->size() == 0 ) ) { // original direction container was not set return; } int c = 0; for(GradientDirectionsContainerType::ConstIterator gdcit = originalDirections->Begin(); gdcit != originalDirections->End(); ++gdcit) { vnl_vector vec = gdcit.Value(); directions->InsertElement(c, vec); c++; } image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::KEEP_ORIGINAL_DIRECTIONS.c_str(), mitk::BoolProperty::New(true) ); image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( directions ) ); } void mitk::DiffusionPropertyHelper::UpdateBValueMap() { BValueMapType b_ValueMap; if(m_Image->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).IsNull()) { } else { b_ValueMap = static_cast(m_Image->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap(); } if(!b_ValueMap.empty()) { b_ValueMap.clear(); } if( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).IsNotNull() ) { GradientDirectionsContainerType::Pointer directions = static_cast( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); GradientDirectionsContainerType::ConstIterator gdcit; for( gdcit = directions->Begin(); gdcit != directions->End(); ++gdcit) { b_ValueMap[GetB_Value(gdcit.Index())].push_back(gdcit.Index()); } } m_Image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str(), mitk::BValueMapProperty::New( b_ValueMap ) ); } bool mitk::DiffusionPropertyHelper::AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision) { GradientDirectionType diff = g1 - g2; GradientDirectionType diff2 = g1 + g2; return diff.two_norm() < precision || diff2.two_norm() < precision; } float mitk::DiffusionPropertyHelper::GetB_Value(unsigned int i) { GradientDirectionsContainerType::Pointer directions = static_cast( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); float b_value = static_cast(m_Image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue(); if(i > directions->Size()-1) return -1; if(directions->ElementAt(i).one_norm() <= 0.0) { return 0; } else { double twonorm = directions->ElementAt(i).two_norm(); double bval = b_value*twonorm*twonorm; if (bval<0) bval = ceil(bval - 0.5); else bval = floor(bval + 0.5); return bval; } } void mitk::DiffusionPropertyHelper::CopyProperties(mitk::Image* source, mitk::Image* target, bool ignore_original_gradients) { mitk::PropertyList::Pointer props = source->GetPropertyList()->Clone(); if (ignore_original_gradients) props->RemoveProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME); target->SetPropertyList(props); } void mitk::DiffusionPropertyHelper::InitializeImage() { if ( m_Image->GetProperty(mitk::DiffusionPropertyHelper::KEEP_ORIGINAL_DIRECTIONS.c_str()).IsNull() ) m_Image->SetProperty( mitk::DiffusionPropertyHelper::KEEP_ORIGINAL_DIRECTIONS.c_str(), mitk::BoolProperty::New(false) ); if( m_Image->GetProperty(mitk::DiffusionPropertyHelper::ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).IsNull() ) { // we don't have the original gradient directions. Therefore use the modified directions and roatate them back. this->UnApplyMeasurementFrameAndRotationMatrix(); } else this->ApplyMeasurementFrameAndRotationMatrix(); this->UpdateBValueMap(); // initialize missing properties mitk::MeasurementFrameProperty::Pointer mf = dynamic_cast( m_Image->GetProperty(MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer()); if( mf.IsNull() ) { //no measurement frame present, identity is assumed MeasurementFrameType identity; identity.set_identity(); m_Image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( identity )); } } bool mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(const mitk::DataNode* node) { if ( node==nullptr ) return false; if ( node->GetData()==nullptr ) return false; return IsDiffusionWeightedImage(dynamic_cast(node->GetData())); } bool mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(const mitk::Image * image) { bool isDiffusionWeightedImage( true ); if( image == nullptr ) { isDiffusionWeightedImage = false; } if( isDiffusionWeightedImage ) { mitk::FloatProperty::Pointer referenceBValue = dynamic_cast(image->GetProperty(REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer()); if( referenceBValue.IsNull() ) { isDiffusionWeightedImage = false; } } unsigned int gradientDirections( 0 ); if( isDiffusionWeightedImage ) { mitk::GradientDirectionsProperty::Pointer gradientDirectionsProperty = dynamic_cast(image->GetProperty(GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer()); if( gradientDirectionsProperty.IsNull() ) { isDiffusionWeightedImage = false; } else { gradientDirections = gradientDirectionsProperty->GetGradientDirectionsContainer()->size(); } } if( isDiffusionWeightedImage ) { unsigned int components = image->GetPixelType().GetNumberOfComponents(); if( components != gradientDirections ) { isDiffusionWeightedImage = false; } } return isDiffusionWeightedImage; } const mitk::DiffusionPropertyHelper::BValueMapType & mitk::DiffusionPropertyHelper::GetBValueMap(const mitk::Image *image) { return dynamic_cast(image->GetProperty(BVALUEMAPPROPERTYNAME.c_str()).GetPointer())->GetBValueMap(); } float mitk::DiffusionPropertyHelper::GetReferenceBValue(const mitk::Image *image) { return dynamic_cast(image->GetProperty(REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer())->GetValue(); } const mitk::DiffusionPropertyHelper::MeasurementFrameType & mitk::DiffusionPropertyHelper::GetMeasurementFrame(const mitk::Image *image) { mitk::MeasurementFrameProperty::Pointer mf = dynamic_cast( image->GetProperty(MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer()); if( mf.IsNull() ) { //no measurement frame present, identity is assumed MeasurementFrameType identity; identity.set_identity(); mf = mitk::MeasurementFrameProperty::New( identity ); } return mf->GetMeasurementFrame(); } mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(const mitk::Image *image) { return dynamic_cast(image->GetProperty(ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer())->GetGradientDirectionsContainer(); } mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::DiffusionPropertyHelper::GetGradientContainer(const mitk::Image *image) { return dynamic_cast(image->GetProperty(GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer())->GetGradientDirectionsContainer(); } bool mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage() const { return IsDiffusionWeightedImage(m_Image); } const mitk::DiffusionPropertyHelper::BValueMapType &mitk::DiffusionPropertyHelper::GetBValueMap() const { return GetBValueMap(m_Image); } float mitk::DiffusionPropertyHelper::GetReferenceBValue() const { return GetReferenceBValue(m_Image); } const mitk::DiffusionPropertyHelper::MeasurementFrameType & mitk::DiffusionPropertyHelper::GetMeasurementFrame() const { return GetMeasurementFrame(m_Image); } mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::DiffusionPropertyHelper::GetOriginalGradientContainer() const { return GetOriginalGradientContainer(m_Image); } mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::DiffusionPropertyHelper::GetGradientContainer() const { return GetGradientContainer(m_Image); } diff --git a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp index 86eadacd24..81da34bea2 100644 --- a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp +++ b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp @@ -1,282 +1,284 @@ /*=================================================================== 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 "mitkFiberBundleVtkReader.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkDiffusionIOMimeTypes.h" mitk::FiberBundleVtkReader::FiberBundleVtkReader() : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_VTK_MIMETYPE_NAME(), "VTK Fiber Bundle Reader" ) { m_ServiceReg = this->RegisterService(); } mitk::FiberBundleVtkReader::FiberBundleVtkReader(const FiberBundleVtkReader &other) :mitk::AbstractFileReader(other) { } mitk::FiberBundleVtkReader * mitk::FiberBundleVtkReader::Clone() const { return new FiberBundleVtkReader(*this); } std::vector > mitk::FiberBundleVtkReader::Read() { std::vector > result; const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::string filename = this->GetInputLocation(); std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); try { MITK_INFO << "Trying to load fiber file as VTK format."; vtkSmartPointer reader = vtkSmartPointer::New(); reader->SetFileName( this->GetInputLocation().c_str() ); if (reader->IsFilePolyData()) { reader->Update(); if ( reader->GetOutput() != nullptr ) { vtkSmartPointer fiberPolyData = reader->GetOutput(); FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData); vtkSmartPointer weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS")); if (weights!=nullptr) { // float weight=0; - for (int i=0; iGetSize(); i++) + for (int i=0; iGetNumberOfValues(); i++) + { if (weights->GetValue(i)<0.0) { MITK_ERROR << "Fiber weight<0 detected! Setting value to 0."; weights->SetValue(i,0); } + } // if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) // { // MITK_INFO << "Weight: " << weights->GetValue(i); // weight = weights->GetValue(i); // } fiberBundle->SetFiberWeights(weights); } vtkSmartPointer fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS")); if (fiberColors!=nullptr) fiberBundle->SetFiberColors(fiberColors); result.push_back(fiberBundle.GetPointer()); return result; } } else MITK_INFO << "File is not VTK format."; } catch(...) { throw; } try { MITK_INFO << "Trying to load fiber file as VTP format."; vtkSmartPointer reader = vtkSmartPointer::New(); reader->SetFileName( this->GetInputLocation().c_str() ); if ( reader->CanReadFile(this->GetInputLocation().c_str()) ) { reader->Update(); if ( reader->GetOutput() != nullptr ) { vtkSmartPointer fiberPolyData = reader->GetOutput(); FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData); vtkSmartPointer weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS")); if (weights!=nullptr) { // float weight=0; // for (int i=0; iGetSize(); i++) // if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) // { // MITK_INFO << "Weight: " << weights->GetValue(i); // weight = weights->GetValue(i); // } fiberBundle->SetFiberWeights(weights); } vtkSmartPointer fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS")); if (fiberColors!=nullptr) fiberBundle->SetFiberColors(fiberColors); result.push_back(fiberBundle.GetPointer()); return result; } } else MITK_INFO << "File is not VTP format."; } catch(...) { throw; } try { MITK_INFO << "Trying to load fiber file as deprecated XML format."; vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); vtkSmartPointer cellArray = vtkSmartPointer::New(); vtkSmartPointer points = vtkSmartPointer::New(); TiXmlDocument doc( this->GetInputLocation().c_str() ); if(doc.LoadFile()) { TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); // save this for later hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("geometry").Element(); // read geometry mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); // read origin mitk::Point3D origin; double temp = 0; pElem->Attribute("origin_x", &temp); origin[0] = temp; pElem->Attribute("origin_y", &temp); origin[1] = temp; pElem->Attribute("origin_z", &temp); origin[2] = temp; geometry->SetOrigin(origin); // read spacing ScalarType spacing[3]; pElem->Attribute("spacing_x", &temp); spacing[0] = temp; pElem->Attribute("spacing_y", &temp); spacing[1] = temp; pElem->Attribute("spacing_z", &temp); spacing[2] = temp; geometry->SetSpacing(spacing); // read transform vtkMatrix4x4* m = vtkMatrix4x4::New(); pElem->Attribute("xx", &temp); m->SetElement(0,0,temp); pElem->Attribute("xy", &temp); m->SetElement(1,0,temp); pElem->Attribute("xz", &temp); m->SetElement(2,0,temp); pElem->Attribute("yx", &temp); m->SetElement(0,1,temp); pElem->Attribute("yy", &temp); m->SetElement(1,1,temp); pElem->Attribute("yz", &temp); m->SetElement(2,1,temp); pElem->Attribute("zx", &temp); m->SetElement(0,2,temp); pElem->Attribute("zy", &temp); m->SetElement(1,2,temp); pElem->Attribute("zz", &temp); m->SetElement(2,2,temp); m->SetElement(0,3,origin[0]); m->SetElement(1,3,origin[1]); m->SetElement(2,3,origin[2]); m->SetElement(3,3,1); geometry->SetIndexToWorldTransformByVtkMatrix(m); // read bounds float bounds[] = {0, 0, 0, 0, 0, 0}; pElem->Attribute("size_x", &temp); bounds[1] = temp; pElem->Attribute("size_y", &temp); bounds[3] = temp; pElem->Attribute("size_z", &temp); bounds[5] = temp; geometry->SetFloatBounds(bounds); geometry->SetImageGeometry(true); pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element(); for( ; pElem ; pElem=pElem->NextSiblingElement()) { TiXmlElement* pElem2 = pElem->FirstChildElement(); vtkSmartPointer container = vtkSmartPointer::New(); for( ; pElem2; pElem2=pElem2->NextSiblingElement()) { Point3D point; pElem2->Attribute("pos_x", &temp); point[0] = temp; pElem2->Attribute("pos_y", &temp); point[1] = temp; pElem2->Attribute("pos_z", &temp); point[2] = temp; geometry->IndexToWorld(point, point); vtkIdType id = points->InsertNextPoint(point.GetDataPointer()); container->GetPointIds()->InsertNextId(id); } cellArray->InsertNextCell(container); } fiberPolyData->SetPoints(points); fiberPolyData->SetLines(cellArray); vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInputData(fiberPolyData); cleaner->Update(); fiberPolyData = cleaner->GetOutput(); FiberBundle::Pointer image = FiberBundle::New(fiberPolyData); result.push_back(image.GetPointer()); return result; } else { MITK_INFO << "File is not deprectaed XML format."; } setlocale(LC_ALL, currLocale.c_str()); MITK_INFO << "Fiber bundle read"; } catch(...) { throw; } throw "Selected file is no vtk readable fiber format (binary or ascii vtk or vtp file)."; return result; } diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h index c3e6c915fc..61e58341b0 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkAstroStickModel.h @@ -1,130 +1,129 @@ /*=================================================================== 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 _MITK_AstroStickModel_H #define _MITK_AstroStickModel_H #include #include namespace mitk { /** * \brief Generates the diffusion signal using a collection of idealised cylinder with zero radius: e^(-bd(ng)²) * */ template< class ScalarType = double > class AstroStickModel : public DiffusionSignalModel< ScalarType > { public: AstroStickModel(); template< class OtherType >AstroStickModel(AstroStickModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_Diffusivity = model->GetDiffusivity(); this->m_Sticks = model->GetSticks(); this->m_NumSticks = model->GetNumSticks(); this->m_RandomizeSticks = model->GetRandomizeSticks(); } ~AstroStickModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ PixelType SimulateMeasurement(); ScalarType SimulateMeasurement(unsigned int dir); void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } - void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } void SetRandomizeSticks(bool randomize=true){ m_RandomizeSticks=randomize; } ///< Random stick configuration in each voxel bool GetRandomizeSticks() { return m_RandomizeSticks; } void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant double GetDiffusivity() { return m_Diffusivity; } void SetNumSticks(unsigned int order) { vnl_matrix sticks; switch (order) { case 1: m_NumSticks = 12; sticks = itk::PointShell<12, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 2: m_NumSticks = 42; sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 3: m_NumSticks = 92; sticks = itk::PointShell<92, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 4: m_NumSticks = 162; sticks = itk::PointShell<162, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; case 5: m_NumSticks = 252; sticks = itk::PointShell<252, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; default: m_NumSticks = 42; sticks = itk::PointShell<42, vnl_matrix_fixed >::DistributePointShell()->as_matrix(); break; } for (int i=0; i namespace mitk { /** * \brief Generates direction independent diffusion measurement employing a scalar diffusion constant d: e^(-bd) * */ template< class ScalarType = double > class BallModel : public DiffusionSignalModel< ScalarType > { public: BallModel(); template< class OtherType >BallModel(BallModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_Diffusivity = model->GetDiffusivity(); } ~BallModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ PixelType SimulateMeasurement(); ScalarType SimulateMeasurement(unsigned int dir); void SetDiffusivity(double D) { m_Diffusivity = D; } double GetDiffusivity() { return m_Diffusivity; } void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } - void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } protected: double m_Diffusivity; ///< Scalar diffusion constant double m_BValue; ///< b-value used to generate the artificial signal }; } #include "mitkBallModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h index 7ff59dfb22..7ad90a30dd 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDiffusionSignalModel.h @@ -1,96 +1,113 @@ /*=================================================================== 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 _MITK_DiffusionSignalModel_H #define _MITK_DiffusionSignalModel_H #include #include #include #include #include #include +#include namespace mitk { /** * \brief Abstract class for diffusion signal models * */ template< class ScalarType = double > class DiffusionSignalModel { public: DiffusionSignalModel() : m_T2(100) , m_T1(0) {} ~DiffusionSignalModel(){} typedef itk::Image ItkDoubleImgType; typedef itk::VariableLengthVector< ScalarType > PixelType; typedef itk::Vector GradientType; typedef std::vector GradientListType; typedef itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRandGenType; + typedef mitk::DiffusionPropertyHelper DPH; /** Realizes actual signal generation. Has to be implemented in subclass. **/ virtual PixelType SimulateMeasurement() = 0; virtual ScalarType SimulateMeasurement(unsigned int dir) = 0; + virtual void SetFiberDirection(GradientType fiberDirection) = 0; GradientType GetFiberDirection(){ return m_FiberDirection; } - virtual void SetGradientList(GradientListType gradientList) = 0; + void SetGradientList(DPH::GradientDirectionsContainerType* gradients) + { + m_GradientList.clear(); + for ( unsigned int i=0; iSize(); ++i ) + { + DPH::GradientDirectionType g_vnl = gradients->GetElement(i); + GradientType g_itk; + g_itk[0] = g_vnl[0]; + g_itk[1] = g_vnl[1]; + g_itk[2] = g_vnl[2]; + m_GradientList.push_back(g_itk); + } + } + + void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } GradientListType GetGradientList(){ return m_GradientList; } GradientType GetGradientDirection(int i) { return m_GradientList.at(i); } void SetT2(double T2) { m_T2 = T2; } double GetT2() { return m_T2; } void SetT1(double T1) { m_T1 = T1; } double GetT1() { return m_T1; } void SetVolumeFractionImage(ItkDoubleImgType::Pointer img){ m_VolumeFractionImage = img; } ItkDoubleImgType::Pointer GetVolumeFractionImage(){ return m_VolumeFractionImage; } void SetRandomGenerator(ItkRandGenType::Pointer randgen){ m_RandGen = randgen; } ItkRandGenType::Pointer GetRandomGenerator(){ return m_RandGen; } void SetSeed(int s) ///< set seed for random generator { if (m_RandGen.IsNull()) m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(s); } unsigned int m_CompartmentId; ///< GUI flag. Which compartment is this model assigned to? protected: GradientType m_FiberDirection; ///< Needed to generate anisotropc signal to determin direction of anisotropy GradientListType m_GradientList; ///< Diffusion gradient direction container double m_T2; ///< Tissue specific transversal relaxation time double m_T1; ///< Tissue specific longitudinal relaxation time ItkDoubleImgType::Pointer m_VolumeFractionImage; ///< Tissue specific volume fraction for each voxel (only relevant for non fiber compartments) ItkRandGenType::Pointer m_RandGen; ///< Random number generator }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h index 382c2dd903..1410833182 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkDotModel.h @@ -1,66 +1,65 @@ /*=================================================================== 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 _MITK_DotModel_H #define _MITK_DotModel_H #include namespace mitk { /** * \brief Generates constant direction independent signal. * */ template< class ScalarType = double > class DotModel : public DiffusionSignalModel< ScalarType > { public: DotModel(); template< class OtherType >DotModel(DotModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); } ~DotModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ PixelType SimulateMeasurement(); ScalarType SimulateMeasurement(unsigned int dir); void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } - void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } protected: }; } #include "mitkDotModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h index 4e6842e05b..3997a72f0a 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkRawShModel.h @@ -1,114 +1,113 @@ /*=================================================================== 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 _MITK_RawShModel_H #define _MITK_RawShModel_H #include #include #include namespace mitk { /** * \brief The spherical harmonic representation of a prototype diffusion weighted MR signal is used to obtain the direction dependent signal. * */ template< class ScalarType = double > class RawShModel : public DiffusionSignalModel< ScalarType > { public: RawShModel(); template< class OtherType >RawShModel(RawShModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_AdcRange = model->GetAdcRange(); this->m_FaRange = model->GetFaRange(); this->m_ShCoefficients = model->GetShCoefficients(); this->m_B0Signal = model->GetB0Signals(); this->m_SphCoords = model->GetSphericalCoordinates(); this->m_ShOrder = model->GetShOrder(); this->m_ModelIndex = model->GetModelIndex(); this->m_MaxNumKernels = model->GetMaxNumKernels(); } ~RawShModel(); typedef itk::Image< double, 3 > ItkDoubleImageType; typedef itk::Image< unsigned char, 3 > ItkUcharImageType; typedef itk::Image< itk::DiffusionTensor3D< double >, 3 > TensorImageType; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter QballFilterType; typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; typedef itk::Matrix< double, 3, 3 > MatrixType; /** Actual signal generation **/ PixelType SimulateMeasurement(); ScalarType SimulateMeasurement(unsigned int dir); bool SetShCoefficients(vnl_vector< double > shCoefficients, double b0); void SetFiberDirection(GradientType fiberDirection); - void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } void SetFaRange(double min, double max){ m_FaRange.first = min; m_FaRange.second = max; } void SetAdcRange(double min, double max){ m_AdcRange.first = min; m_AdcRange.second = max; } void SetMaxNumKernels(unsigned int max){ m_MaxNumKernels = max; } unsigned int GetNumberOfKernels(); std::pair< double, double > GetFaRange(){ return m_FaRange; } std::pair< double, double > GetAdcRange(){ return m_AdcRange; } unsigned int GetMaxNumKernels(){ return m_MaxNumKernels; } void Clear(); std::vector< vnl_vector< double > > GetShCoefficients(){ return m_ShCoefficients; } std::vector< double > GetB0Signals(){ return m_B0Signal; } vnl_matrix GetSphericalCoordinates(){ return m_SphCoords; } unsigned int GetShOrder(){ return m_ShOrder; } int GetModelIndex(){ return m_ModelIndex; } double GetBaselineSignal(int index){ return m_B0Signal.at(index); } vnl_vector< double > GetCoefficients(int listIndex){ return m_ShCoefficients.at(listIndex); } bool SampleKernels(Image::Pointer diffImg, ItkUcharImageType::Pointer maskImage, TensorImageType::Pointer tensorImage=nullptr, QballFilterType::CoefficientImageType::Pointer coeffImage=nullptr, ItkDoubleImageType::Pointer adcImage=nullptr); protected: void Cart2Sph( GradientListType gradients ); void RandomModel(); std::vector< vnl_vector< double > > m_ShCoefficients; std::vector< double > m_B0Signal; std::vector< GradientType > m_PrototypeMaxDirection; vnl_matrix m_SphCoords; std::pair< double, double > m_AdcRange; std::pair< double, double > m_FaRange; unsigned int m_ShOrder; int m_ModelIndex; unsigned int m_MaxNumKernels; }; } #include "mitkRawShModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h index 7ecf8635d7..bc61237f84 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkStickModel.h @@ -1,76 +1,75 @@ /*=================================================================== 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 _MITK_StickModel_H #define _MITK_StickModel_H #include namespace mitk { /** * \brief Generates the diffusion signal using an idealised cylinder with zero radius: e^(-bd(ng)²) * */ template< class ScalarType = double > class StickModel : public DiffusionSignalModel< ScalarType > { public: StickModel(); template< class OtherType >StickModel(StickModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_Diffusivity = model->GetDiffusivity(); } ~StickModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ PixelType SimulateMeasurement(); ScalarType SimulateMeasurement(unsigned int dir); void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } void SetDiffusivity(double diffusivity) { m_Diffusivity = diffusivity; } ///< Scalar diffusion constant double GetDiffusivity() { return m_Diffusivity; } void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } - void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } protected: double m_Diffusivity; ///< Scalar diffusion constant double m_BValue; ///< b-value used to generate the artificial signal }; } #include "mitkStickModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h index 7734ba3cc3..4d9e73a953 100644 --- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.h +++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkTensorModel.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 _MITK_TensorModel_H #define _MITK_TensorModel_H #include #include namespace mitk { /** * \brief Generates diffusion measurement employing a second rank tensor model: e^(-bg^TDg) * */ template< class ScalarType = double > class TensorModel : public DiffusionSignalModel< ScalarType > { public: TensorModel(); template< class OtherType >TensorModel(TensorModel* model) { this->m_CompartmentId = model->m_CompartmentId; this->m_T2 = model->GetT2(); this->m_FiberDirection = model->GetFiberDirection(); this->m_GradientList = model->GetGradientList(); this->m_VolumeFractionImage = model->GetVolumeFractionImage(); this->m_RandGen = model->GetRandomGenerator(); this->m_BValue = model->GetBvalue(); this->m_KernelDirection = model->GetKernelDirection(); this->m_KernelTensorMatrix = model->GetKernelTensorMatrix(); } ~TensorModel(); typedef typename DiffusionSignalModel< ScalarType >::PixelType PixelType; typedef itk::DiffusionTensor3D< ScalarType > ItkTensorType; typedef typename DiffusionSignalModel< ScalarType >::GradientType GradientType; typedef typename DiffusionSignalModel< ScalarType >::GradientListType GradientListType; /** Actual signal generation **/ PixelType SimulateMeasurement(); ScalarType SimulateMeasurement(unsigned int dir); void SetBvalue(double bValue) { m_BValue = bValue; } ///< b-value used to generate the artificial signal double GetBvalue() { return m_BValue; } void SetDiffusivity1(double d1){ m_KernelTensorMatrix[0][0] = d1; } void SetDiffusivity2(double d2){ m_KernelTensorMatrix[1][1] = d2; } void SetDiffusivity3(double d3){ m_KernelTensorMatrix[2][2] = d3; } double GetDiffusivity1() { return m_KernelTensorMatrix[0][0]; } double GetDiffusivity2() { return m_KernelTensorMatrix[1][1]; } double GetDiffusivity3() { return m_KernelTensorMatrix[2][2]; } void SetFiberDirection(GradientType fiberDirection){ this->m_FiberDirection = fiberDirection; } - void SetGradientList(GradientListType gradientList) { this->m_GradientList = gradientList; } GradientType GetKernelDirection(){ return m_KernelDirection; } vnl_matrix_fixed GetKernelTensorMatrix(){ return m_KernelTensorMatrix; } protected: /** Calculates tensor matrix from FA and ADC **/ void UpdateKernelTensor(); GradientType m_KernelDirection; ///< Direction of the kernel tensors principal eigenvector vnl_matrix_fixed m_KernelTensorMatrix; ///< 3x3 matrix containing the kernel tensor values double m_BValue; ///< b-value used to generate the artificial signal }; } #include "mitkTensorModel.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp index 159b51cc8a..d53e7b4d7e 100755 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp @@ -1,2325 +1,2325 @@ /*=================================================================== 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 _USE_MATH_DEFINES #include "mitkFiberBundle.h" #include #include #include #include "mitkImagePixelReadAccessor.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs"; using namespace std; mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData ) : m_NumFibers(0) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != nullptr) m_FiberPolyData = fiberPolyData; this->UpdateFiberGeometry(); this->GenerateFiberIds(); this->ColorFibersByOrientation(); } mitk::FiberBundle::~FiberBundle() { } mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy() { mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); newFib->SetFiberColors(this->m_FiberColors); newFib->SetFiberWeights(this->m_FiberWeights); return newFib; } vtkSmartPointer mitk::FiberBundle::GeneratePolyDataByIds(std::vector fiberIds) { vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); auto finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*finIt < 0 || *finIt>GetNumFibers()){ MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; break; } vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); vtkSmartPointer fibPoints = fiber->GetPoints(); vtkSmartPointer newFiber = vtkSmartPointer::New(); newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); for(int i=0; iGetNumberOfPoints(); i++) { newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); } newLineSet->InsertNextCell(newFiber); ++finIt; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib) { if (fib==nullptr) { MITK_WARN << "trying to call AddBundle with nullptr argument"; return nullptr; } MITK_INFO << "Adding fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); unsigned int counter = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // add new fiber bundle for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, fib->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // subtract two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib) { MITK_INFO << "Subtracting fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector< std::vector< itk::Point > > points1; for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = GetItkPoint(points->GetPoint(0)); itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); points1.push_back( {start, end} ); } std::vector< std::vector< itk::Point > > points2; for( int i=0; iGetNumFibers(); i++ ) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = GetItkPoint(points->GetPoint(0)); itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); points2.push_back( {start, end} ); } int progress = 0; std::vector< int > ids; #pragma omp parallel for for (int i=0; i<(int)points1.size(); i++) { #pragma omp critical { progress++; std::cout << (int)(100*(float)progress/points1.size()) << "%" << '\r'; cout.flush(); } bool match = false; for (unsigned int j=0; jGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(points->GetPoint(j)); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if(vNewLines->GetNumberOfCells()==0) return nullptr; // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle return mitk::FiberBundle::New(vNewPolyData); } itk::Point mitk::FiberBundle::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } /* * set PolyData (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { if (fiberPD == nullptr) this->m_FiberPolyData = vtkSmartPointer::New(); else m_FiberPolyData->DeepCopy(fiberPD); m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (updateGeometry) UpdateFiberGeometry(); GenerateFiberIds(); ColorFibersByOrientation(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundle::GetFiberPolyData() const { return m_FiberPolyData; } void mitk::FiberBundle::ColorFibersByOrientation() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also PolyData needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= vtkPoints* extrPoints = nullptr; extrPoints = m_FiberPolyData->GetPoints(); int numOfPoints = 0; if (extrPoints!=nullptr) numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = 4; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(numOfPoints * componentSize); m_FiberColors->SetNumberOfComponents(componentSize); m_FiberColors->SetName("FIBER_COLORS"); int numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) return; /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); /* single fiber checkpoints: is number of points valid */ if (pointsPerFiber > 1) { /* operate on points of single fiber */ for (int i=0; i 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); rgba[3] = (unsigned char) (255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); rgba[3] = (unsigned char) (255.0); } m_FiberColors->InsertTypedTuple(idList[i], rgba); } } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; } else { MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByCurvature(bool, bool normalize) { double window = 5; //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = 4; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize); m_FiberColors->SetNumberOfComponents(componentSize); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); vector< double > values; double min = 1; double max = 0; MITK_INFO << "Coloring fibers by curvature"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures for (int j=0; j > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0); while(dist1) { double p1[3]; points->GetPoint(c-1, p1); double p2[3]; points->GetPoint(c, p2); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==j) meanV += v; c--; } c = j; dist = 0; while(distGetPoint(c, p1); double p2[3]; points->GetPoint(c+1, p2); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==j) meanV += v; c++; } meanV.normalize(); double dev = 0; for (unsigned int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); dev = 1.0-dev/180.0; values.push_back(dev); if (devmax) max = dev; } } unsigned int count = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); for (int j=0; j1) dev = 1; lookupTable->GetColor(dev, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba); count++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray) { for(long i=0; iGetNumberOfTuples(); i++) { double faValue = FAValArray->GetValue(i); faValue = faValue * 255.0; m_FiberColors->SetComponent(i,3, (unsigned char) faValue ); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ResetFiberOpacity() { for(long i=0; iGetNumberOfTuples(); i++) m_FiberColors->SetComponent(i,3, 255.0 ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize) { mitkPixelTypeMultiplex3( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } template void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); unsigned char rgba[4] = {0,0,0,0}; vtkPoints* pointSet = m_FiberPolyData->GetPoints(); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); double min = 9999999; double max = -9999999; for(long i=0; iGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; double pixelValue = readimage.GetPixelByWorldCoordinates(px); if (pixelValue>max) max = pixelValue; if (pixelValueGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; double pixelValue = readimage.GetPixelByWorldCoordinates(px); if (normalize) pixelValue = (pixelValue-min)/(max-min); else if (pixelValue>1) pixelValue = 1; double color[3]; lookupTable->GetColor(1-pixelValue, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); if (opacity) rgba[3] = (unsigned char) (255.0 * pixelValue); else rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByFiberWeights(bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); unsigned char rgba[4] = {0,0,0,0}; unsigned int counter = 0; float max = -999999; float min = 999999; for (int i=0; iGetFiberWeight(i); if (weight>max) max = weight; if (weightGetCell(i); int numPoints = cell->GetNumberOfPoints(); double weight = this->GetFiberWeight(i); for (int j=0; j1) v = 1; double color[3]; lookupTable->GetColor(1-v, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); if (opacity) rgba[3] = (unsigned char) (255.0 * v); else rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTypedTuple(counter, rgba); counter++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); unsigned char rgba[4] = {0,0,0,0}; for(long i=0; iGetNumberOfPoints(); ++i) { rgba[0] = (unsigned char) r; rgba[1] = (unsigned char) g; rgba[2] = (unsigned char) b; rgba[3] = (unsigned char) alpha; m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::GenerateFiberIds() { if (m_FiberPolyData == nullptr) return; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInputData(m_FiberPolyData); idFiberFilter->CellIdsOn(); // idFiberFilter->PointIdsOn(); // point id's are not needed idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); idFiberFilter->FieldDataOn(); idFiberFilter->Update(); m_FiberIdDataSet = idFiberFilter->GetOutput(); } -mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds, float fraction) +mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds, float fraction, bool do_resampling) { vtkSmartPointer PolyData = m_FiberPolyData; - if (anyPoint) + if (anyPoint && do_resampling) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); fibCopy->ResampleLinear(minSpacing/5); PolyData = fibCopy->GetFiberPolyData(); } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Extracting fibers with mask image"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cellOriginal = m_FiberPolyData->GetCell(i); int numPointsOriginal = cellOriginal->GetNumberOfPoints(); vtkPoints* pointsOriginal = cellOriginal->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); if (numPoints>1 && numPointsOriginal) { if (anyPoint) { int inside = 0; int outside = 0; if (!invert) { for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 ) { inside++; if (fraction==0) break; } else outside++; } float current_fraction = 0.0; if (inside+outside>0) current_fraction = (float)inside/(inside+outside); if (current_fraction>fraction) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else { bool includeFiber = true; for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) ) { inside++; includeFiber = false; break; } else outside++; } if (includeFiber) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } else { double* start = pointsOriginal->GetPoint(0); itk::Point itkStart; itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; itk::Index<3> idxStart; mask->TransformPhysicalPointToIndex(itkStart, idxStart); double* end = pointsOriginal->GetPoint(numPointsOriginal-1); itk::Point itkEnd; itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; itk::Index<3> idxEnd; mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); if (invert) { if (bothEnds) { if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else { if (bothEnds) { if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } } vtkNewCells->InsertNextCell(container); } if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); return mitk::FiberBundle::New(newPolyData); } mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); fibCopy->ResampleLinear(minSpacing/10); vtkSmartPointer PolyData =fibCopy->GetFiberPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Cutting fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); if (numPoints>1) { int newNumPoints = 0; for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if (newNumPoints>0) { vtkNewCells->InsertNextCell(container); newNumPoints = 0; container = vtkSmartPointer::New(); } } if (newNumPoints>0) vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData); newFib->Compress(0.1); return newFib; } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage) { if (roi==nullptr || !(dynamic_cast(roi->GetData()) || dynamic_cast(roi->GetData())) ) return nullptr; std::vector tmp = ExtractFiberIdSubset(roi, storage); if (tmp.size()<=0) return mitk::FiberBundle::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); return mitk::FiberBundle::New(pTmp); } std::vector mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage) { std::vector result; if (roi==nullptr || roi->GetData()==nullptr) return result; mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); if (!pfc.IsNull()) // handle composite { DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); if (children->size()==0) return result; switch (pfc->getOperationType()) { case 0: // AND { MITK_INFO << "AND"; result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(std::min(result.size(),inRoi.size())); it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( it - rest.begin() ); result = rest; } break; } case 1: // OR { MITK_INFO << "OR"; result = ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { it = result.end(); std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); result.insert(it, inRoi.begin(), inRoi.end()); } // remove duplicates sort(result.begin(), result.end()); it = unique(result.begin(), result.end()); result.resize( it - result.begin() ); break; } case 2: // NOT { MITK_INFO << "NOT"; for(long i=0; iGetNumFibers(); i++) result.push_back(i); std::vector::iterator it; for (unsigned int i=0; iSize(); ++i) { std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(result.size()-inRoi.size()); it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( it - rest.begin() ); result = rest; } break; } } } else if ( dynamic_cast(roi->GetData()) ) // actual extraction { if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarPoly = dynamic_cast(roi->GetData()); //create vtkPolygon using controlpoints from planarFigure polygon vtkSmartPointer polygonVtk = vtkSmartPointer::New(); for (unsigned int i=0; iGetNumberOfControlPoints(); ++i) { itk::Point p = planarPoly->GetWorldControlPoint(i); vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] ); polygonVtk->GetPointIds()->InsertNextId(id); } MITK_INFO << "Extracting with polygon"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); double tolerance = 0.001; // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection double pcoords[3] = {0,0,0}; int subId = 0; int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId); if (iD!=0) { result.push_back(i); break; } } } } else if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(roi->GetData()); Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal(); planeNormal.Normalize(); //calculate circle radius mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint double radius = V1w.EuclideanDistanceTo(V2w); radius *= radius; MITK_INFO << "Extracting with circle"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x); if (iD!=0) { double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]); if( dist <= radius) { result.push_back(i); break; } } } } } return result; } return result; } void mitk::FiberBundle::UpdateFiberGeometry() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInputData(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = m_FiberPolyData->GetNumberOfCells(); if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints()) this->ColorFibersByOrientation(); - if (m_FiberWeights->GetSize()!=m_NumFibers) + if (m_FiberWeights->GetNumberOfValues()!=m_NumFibers) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberWeights->SetNumberOfValues(m_NumFibers); this->SetFiberWeights(1); } if (m_NumFibers<=0) // no fibers present; apply default geometry { m_MinFiberLength = 0; m_MaxFiberLength = 0; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(false); float b[] = {0, 1, 0, 1, 0, 1}; geometry->SetFloatBounds(b); SetGeometry(geometry); return; } double b[6]; m_FiberPolyData->GetBounds(b); // calculate statistics for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); length += dist; } m_FiberLengths.push_back(length); m_MeanFiberLength += length; if (i==0) { m_MinFiberLength = length; m_MaxFiberLength = length; } else { if (lengthm_MaxFiberLength) m_MaxFiberLength = length; } } m_MeanFiberLength /= m_NumFibers; std::vector< float > sortedLengths = m_FiberLengths; std::sort(sortedLengths.begin(), sortedLengths.end()); for (int i=0; i1) m_LengthStDev /= (m_NumFibers-1); else m_LengthStDev = 0; m_LengthStDev = std::sqrt(m_LengthStDev); m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetFloatBounds(b); this->SetGeometry(geometry); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const { return m_FiberWeights->GetValue(fiber); } void mitk::FiberBundle::SetFiberWeights(float newWeight) { - for (int i=0; iGetSize(); i++) + for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, newWeight); } void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer weights) { - if (m_NumFibers!=weights->GetSize()) + if (m_NumFibers!=weights->GetNumberOfValues()) { - MITK_INFO << "Weights array not equal to number of fibers!"; + MITK_INFO << "Weights array not equal to number of fibers! " << weights->GetNumberOfValues() << " vs " << m_NumFibers; return; } - for (int i=0; iGetSize(); i++) + for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, weights->GetValue(i)); m_FiberWeights->SetName("FIBER_WEIGHTS"); } void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight) { m_FiberWeights->SetValue(fiber, weight); } void mitk::FiberBundle::SetFiberColors(vtkSmartPointer fiberColors) { for(long i=0; iGetNumberOfPoints(); ++i) { unsigned char source[4] = {0,0,0,0}; fiberColors->GetTypedTuple(i, source); unsigned char target[4] = {0,0,0,0}; target[0] = source[0]; target[1] = source[1]; target[2] = source[2]; target[3] = source[3]; m_FiberColors->InsertTypedTuple(i, target); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; m = rot*m; return m; } itk::Point mitk::FiberBundle::TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); point[0] -= center[0]; point[1] -= center[1]; point[2] -= center[2]; point = rot*point; point[0] += center[0]+tx; point[1] += center[1]+ty; point[2] += center[2]+tz; itk::Point out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2]; return out; } void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rot*dir; dir[0] += center[0]+tx; dir[1] += center[1]+ty; dir[2] += center[2]+tz; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z) { x = x*M_PI/180; y = y*M_PI/180; z = z*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rotZ*rotY*rotX*dir; dir[0] += center[0]; dir[1] += center[1]; dir[2] += center[2]; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter) { MITK_INFO << "Scaling fibers"; boost::progress_display disp(m_NumFibers); mitk::BaseGeometry* geom = this->GetGeometry(); mitk::Point3D c = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); if (subtractCenter) { p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; } p[0] *= x; p[1] *= y; p[2] *= z; if (subtractCenter) { p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; } vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TranslateFibers(double x, double y, double z) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[0] += x; p[1] += y; p[2] += z; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::MirrorFibers(unsigned int axis) { if (axis>2) return; MITK_INFO << "Mirroring fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[axis] = -p[axis]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RemoveDir(vnl_vector_fixed dir, double threshold) { dir.normalize(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); bool discard = false; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); vnl_vector_fixed< double, 3 > v1; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; if (v1.magnitude()>0.001) { v1.normalize(); if (fabs(dot_product(v1,dir))>threshold) { discard = true; break; } } } if (!discard) { for (int j=0; jGetPoint(j, p1); vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); // UpdateColorCoding(); // UpdateFiberGeometry(); } bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double p3[3]; points->GetPoint(j+2, p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; v2[0] = p3[0]-p2[0]; v2[1] = p3[1]-p2[1]; v2[2] = p3[2]-p2[2]; v3[0] = p1[0]-p3[0]; v3[1] = p1[1]-p3[1]; v3[2] = p1[2]-p3[2]; float a = v1.magnitude(); float b = v2.magnitude(); float c = v3.magnitude(); float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); if (deleteFibers && rInsertNextCell(container); container = vtkSmartPointer::New(); } else if (j==numPoints-3) { id = vtkNewPoints->InsertNextPoint(p2); container->GetPointIds()->InsertNextId(id); id = vtkNewPoints->InsertNextPoint(p3); container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); } } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM) { MITK_INFO << "Removing short fibers"; if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers { MITK_WARN << "Process aborted. No fibers would be left!"; return false; } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); float min = m_MaxFiberLength; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); if (m_FiberLengths.at(i)GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Removing long fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias ) { if (pointDistance<=0) return; vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::New(); //in smoothpoints the interpolated points representing a fiber are stored. //in vtkcells all polylines are stored, actually all id's of them are stored vtkSmartPointer vtkSmoothCells = vtkSmartPointer::New(); //cellcontainer for smoothed lines vtkIdType pointHelperCnt = 0; MITK_INFO << "Smoothing fibers"; vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); boost::progress_display disp(m_NumFibers); #pragma omp parallel for for (int i=0; i newPoints = vtkSmartPointer::New(); float length = 0; float weight = 1; #pragma omp critical { length = m_FiberLengths.at(i); weight = m_FiberWeights->GetValue(i); ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jInsertNextPoint(points->GetPoint(j)); } int sampling = std::ceil(length/pointDistance); vtkSmartPointer xSpline = vtkSmartPointer::New(); vtkSmartPointer ySpline = vtkSmartPointer::New(); vtkSmartPointer zSpline = vtkSmartPointer::New(); xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); vtkSmartPointer spline = vtkSmartPointer::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(newPoints); vtkSmartPointer functionSource = vtkSmartPointer::New(); functionSource->SetParametricFunction(spline); functionSource->SetUResolution(sampling); functionSource->SetVResolution(sampling); functionSource->SetWResolution(sampling); functionSource->Update(); vtkPolyData* outputFunction = functionSource->GetOutput(); vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber vtkSmartPointer smoothLine = vtkSmartPointer::New(); smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints()); #pragma omp critical { for (int j=0; jGetNumberOfPoints(); j++) { smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt); vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); } newFiberWeights->SetValue(vtkSmoothCells->GetNumberOfCells(), weight); vtkSmoothCells->InsertNextCell(smoothLine); pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints(); } } SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ResampleSpline(float pointDistance) { ResampleSpline(pointDistance, 0, 0, 0 ); } unsigned long mitk::FiberBundle::GetNumberOfPoints() const { unsigned long points = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); points += cell->GetNumberOfPoints(); } return points; } void mitk::FiberBundle::Compress(float error) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Compressing fibers"; unsigned long numRemovedPoints = 0; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } // calculate curvatures int numPoints = vertices.size(); std::vector< int > removedPoints; removedPoints.resize(numPoints, 0); removedPoints[0]=-1; removedPoints[numPoints-1]=-1; vtkSmartPointer container = vtkSmartPointer::New(); int remCounter = 0; bool pointFound = true; while (pointFound) { pointFound = false; double minError = error; int removeIndex = -1; for (unsigned int j=0; j candV = vertices.at(j); int validP = -1; vnl_vector_fixed< double, 3 > pred; for (int k=j-1; k>=0; k--) if (removedPoints[k]<=0) { pred = vertices.at(k); validP = k; break; } int validS = -1; vnl_vector_fixed< double, 3 > succ; for (int k=j+1; k=0 && validS>=0) { double a = (candV-pred).magnitude(); double b = (candV-succ).magnitude(); double c = (pred-succ).magnitude(); double s=0.5*(a+b+c); double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c))); if (hcInsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); numRemovedPoints += remCounter; vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { MITK_INFO << "Removed points: " << numRemovedPoints; SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } void mitk::FiberBundle::ResampleLinear(double pointDistance) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Resampling fibers (linear)"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); container->GetPointIds()->InsertNextId(id); } for (unsigned int j=1; j vec = vertices.at(j) - lastV; double new_dist = vec.magnitude(); if (new_dist >= pointDistance) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-pointDistance <= mitk::eps ) { vec.normalize(); newV += vec * pointDistance; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if (j==vertices.size()-1 && new_dist>0.0001) { #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } // reapply selected colorcoding in case PolyData structure has changed bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps) { if (fib==nullptr) { MITK_INFO << "Reference bundle is nullptr!"; return false; } if (m_NumFibers!=fib->GetNumFibers()) { MITK_INFO << "Unequal number of fibers!"; MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); return false; } for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i); int numPoints2 = cell2->GetNumberOfPoints(); vtkPoints* points2 = cell2->GetPoints(); if (numPoints2!=numPoints) { MITK_INFO << "Unequal number of points in fiber " << i << "!"; MITK_INFO << numPoints2 << " vs. " << numPoints; return false; } for (int j=0; jGetPoint(j); double* p2 = points2->GetPoint(j); if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps) { MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!"; MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2]; MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2]; return false; } } } return true; } void mitk::FiberBundle::PrintSelf(std::ostream &os, itk::Indent indent) const { os << indent << this->GetNameOfClass() << ":\n"; os << indent << "Number of fibers: " << this->GetNumFibers() << std::endl; os << indent << "Min. fiber length: " << this->GetMinFiberLength() << std::endl; os << indent << "Max. fiber length: " << this->GetMaxFiberLength() << std::endl; os << indent << "Mean fiber length: " << this->GetMeanFiberLength() << std::endl; os << indent << "Median fiber length: " << this->GetMedianFiberLength() << std::endl; os << indent << "STDEV fiber length: " << this->GetLengthStDev() << std::endl; os << indent << "Number of points: " << this->GetNumberOfPoints() << std::endl; Superclass::PrintSelf(os, indent); } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundle::UpdateOutputInformation() { } void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundle::VerifyRequestedRegion() { return true; } void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* ) { } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h index 10d2264cf5..a49a44e2dd 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h @@ -1,179 +1,179 @@ /*=================================================================== 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 _MITK_FiberBundle_H #define _MITK_FiberBundle_H //includes for MITK datastructure #include #include #include #include #include #include #include //includes storing fiberdata #include #include #include #include #include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class MITKFIBERTRACKING_EXPORT FiberBundle : public BaseData { public: typedef itk::Image ItkUcharImgType; // fiber colorcodings static const char* FIBER_ID_ARRAY; virtual void UpdateOutputInformation() override; virtual void SetRequestedRegionToLargestPossibleRegion() override; virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override; virtual bool VerifyRequestedRegion() override; virtual void SetRequestedRegion(const itk::DataObject*) override; mitkClassMacro( FiberBundle, BaseData ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(Self, vtkSmartPointer) // custom constructor // colorcoding related methods void ColorFibersByFiberWeights(bool opacity, bool normalize); void ColorFibersByCurvature(bool opacity, bool normalize); void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity, bool normalize); template void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity, bool normalize); void ColorFibersByOrientation(); void SetFiberOpacity(vtkDoubleArray *FAValArray); void ResetFiberOpacity(); void SetFiberColors(vtkSmartPointer fiberColors); void SetFiberColors(float r, float g, float b, float alpha=255); vtkSmartPointer GetFiberColors() const { return m_FiberColors; } // fiber compression void Compress(float error = 0.0); // fiber resampling void ResampleSpline(float pointDistance=1); void ResampleSpline(float pointDistance, double tension, double continuity, double bias ); void ResampleLinear(double pointDistance=1); bool RemoveShortFibers(float lengthInMM); bool RemoveLongFibers(float lengthInMM); bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers); void MirrorFibers(unsigned int axis); void RotateAroundAxis(double x, double y, double z); void TranslateFibers(double x, double y, double z); void ScaleFibers(double x, double y, double z, bool subtractCenter=true); void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz); void RemoveDir(vnl_vector_fixed dir, double threshold); itk::Point TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz); itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz); // add/subtract fibers FiberBundle::Pointer AddBundle(FiberBundle* fib); FiberBundle::Pointer SubtractBundle(FiberBundle* fib); // fiber subset extraction FiberBundle::Pointer ExtractFiberSubset(DataNode *roi, DataStorage* storage); std::vector ExtractFiberIdSubset(DataNode* roi, DataStorage* storage); - FiberBundle::Pointer ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert=false, bool bothEnds=true, float fraction=0.0); + FiberBundle::Pointer ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert=false, bool bothEnds=true, float fraction=0.0, bool do_resampling=true); FiberBundle::Pointer RemoveFibersOutside(ItkUcharImgType* mask, bool invert=false); vtkSmartPointer GeneratePolyDataByIds( std::vector ); // TODO: make protected void GenerateFiberIds(); // TODO: make protected // get/set data vtkSmartPointer GetFiberWeights() const { return m_FiberWeights; } float GetFiberWeight(unsigned int fiber) const; void SetFiberWeights(float newWeight); void SetFiberWeight(unsigned int fiber, float weight); void SetFiberWeights(vtkSmartPointer weights); void SetFiberPolyData(vtkSmartPointer, bool updateGeometry = true); vtkSmartPointer GetFiberPolyData() const; itkGetConstMacro( NumFibers, int) //itkGetMacro( FiberSampling, int) itkGetConstMacro( MinFiberLength, float ) itkGetConstMacro( MaxFiberLength, float ) itkGetConstMacro( MeanFiberLength, float ) itkGetConstMacro( MedianFiberLength, float ) itkGetConstMacro( LengthStDev, float ) itkGetConstMacro( UpdateTime2D, itk::TimeStamp ) itkGetConstMacro( UpdateTime3D, itk::TimeStamp ) void RequestUpdate2D(){ m_UpdateTime2D.Modified(); } void RequestUpdate3D(){ m_UpdateTime3D.Modified(); } void RequestUpdate(){ m_UpdateTime2D.Modified(); m_UpdateTime3D.Modified(); } unsigned long GetNumberOfPoints() const; // copy fiber bundle mitk::FiberBundle::Pointer GetDeepCopy(); // compare fiber bundles bool Equals(FiberBundle* fib, double eps=0.01); itkSetMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) itkGetConstMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) protected: FiberBundle( vtkPolyData* fiberPolyData = nullptr ); virtual ~FiberBundle(); itk::Point GetItkPoint(double point[3]); // calculate geometry from fiber extent void UpdateFiberGeometry(); virtual void PrintSelf(std::ostream &os, itk::Indent indent) const override; private: // actual fiber container vtkSmartPointer m_FiberPolyData; // contains fiber ids vtkSmartPointer m_FiberIdDataSet; int m_NumFibers; vtkSmartPointer m_FiberColors; vtkSmartPointer m_FiberWeights; std::vector< float > m_FiberLengths; float m_MinFiberLength; float m_MaxFiberLength; float m_MeanFiberLength; float m_MedianFiberLength; float m_LengthStDev; itk::TimeStamp m_UpdateTime2D; itk::TimeStamp m_UpdateTime3D; mitk::BaseGeometry::Pointer m_ReferenceGeometry; }; } // namespace mitk #endif /* _MITK_FiberBundle_H */ diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/Fiberfox.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/Fiberfox.cpp index 77d6388a9c..30fd887c00 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/Fiberfox.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/Fiberfox.cpp @@ -1,214 +1,209 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include using namespace mitk; /*! * \brief Command line interface to Fiberfox. * Simulate a diffusion-weighted image from a tractogram using the specified parameter file. */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiberfox"); parser.setCategory("Diffusion Simulation Tools"); parser.setContributor("MIC"); - parser.setDescription("Command line interface to Fiberfox." - " Simulate a diffusion-weighted image from a tractogram using the specified parameter file."); + parser.setDescription("Command line interface to Fiberfox." " Simulate a diffusion-weighted image from a tractogram using the specified parameter file."); parser.setArgumentPrefix("--", "-"); - parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output root:", - "output root", us::Any(), false); - parser.addArgument("parameters", "p", mitkCommandLineParser::InputFile, "Parameter file:", - "fiberfox parameter file (.ffp)", us::Any(), false); - parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", - "Input tractogram or diffusion-weighted image.", us::Any(), false); - parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Output additional images:", - "output volume fraction images etc.", us::Any()); + parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output root:", "output root", us::Any(), false); + parser.addArgument("parameters", "p", mitkCommandLineParser::InputFile, "Parameter file:", "fiberfox parameter file (.ffp)", us::Any(), false); + parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "Input tractogram or diffusion-weighted image.", us::Any(), false); + parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Output additional images:", "output volume fraction images etc.", us::Any()); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) { return EXIT_FAILURE; } string outName = us::any_cast(parsedArgs["out"]); string paramName = us::any_cast(parsedArgs["parameters"]); string input=""; if (parsedArgs.count("input")) { input = us::any_cast(parsedArgs["input"]); } bool verbose = false; if (parsedArgs.count("verbose")) { verbose = us::any_cast(parsedArgs["verbose"]); } FiberfoxParameters parameters; parameters.LoadParameters(paramName); // Test if /path/dir is an existing directory: string file_extension = ""; if( itksys::SystemTools::FileIsDirectory( outName ) ) { while( *(--(outName.cend())) == '/') { outName.pop_back(); } outName = outName + '/'; parameters.m_Misc.m_OutputPath = outName; outName = outName + parameters.m_Misc.m_OutputPrefix; // using default m_OutputPrefix as initialized. } else { // outName is NOT an existing directory, so we need to remove all trailing slashes: while( *(--(outName.cend())) == '/') { outName.pop_back(); } // now split up the given outName into directory and (prefix of) filename: if( ! itksys::SystemTools::GetFilenamePath( outName ).empty() && itksys::SystemTools::FileIsDirectory(itksys::SystemTools::GetFilenamePath( outName ) ) ) { parameters.m_Misc.m_OutputPath = itksys::SystemTools::GetFilenamePath( outName ) + '/'; } else { parameters.m_Misc.m_OutputPath = itksys::SystemTools::GetCurrentWorkingDirectory() + '/'; } file_extension = itksys::SystemTools::GetFilenameExtension(outName); if( ! itksys::SystemTools::GetFilenameWithoutExtension( outName ).empty() ) { parameters.m_Misc.m_OutputPrefix = itksys::SystemTools::GetFilenameWithoutExtension( outName ); } else { parameters.m_Misc.m_OutputPrefix = "fiberfox"; } outName = parameters.m_Misc.m_OutputPath + parameters.m_Misc.m_OutputPrefix; } // check if log file already exists and avoid overwriting existing files: std::string NameTest = outName; int c = 0; while( itksys::SystemTools::FileExists( outName + ".log" ) && c <= std::numeric_limits::max() ) { outName = NameTest + "_" + boost::lexical_cast(c); ++c; } if (verbose) { MITK_DEBUG << outName << ".ffp"; parameters.SaveParameters(outName+".ffp"); } mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "Fiberbundles"}, {}); mitk::BaseData::Pointer inputData = mitk::IOUtil::Load(input, &functor)[0]; itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); if ( dynamic_cast(inputData.GetPointer()) ) // simulate dataset from fibers { tractsToDwiFilter->SetFiberBundle(dynamic_cast(inputData.GetPointer())); } else if ( dynamic_cast(inputData.GetPointer()) ) // add artifacts to existing image { typedef itk::VectorImage< short, 3 > ItkDwiType; mitk::Image::Pointer diffImg = dynamic_cast(inputData.GetPointer()); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); parameters.m_SignalGen.m_SignalScale = 1; parameters.m_SignalGen.m_ImageRegion = itkVectorImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkVectorImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkVectorImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkVectorImagePointer->GetDirection(); parameters.m_SignalGen.m_Bvalue = static_cast (diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() ) ->GetValue(); parameters.m_SignalGen.SetGradienDirections( static_cast ( diffImg->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ) ->GetGradientDirectionsContainer() ); tractsToDwiFilter->SetInputImage(itkVectorImagePointer); } tractsToDwiFilter->SetParameters(parameters); tractsToDwiFilter->Update(); mitk::Image::Pointer image = mitk::GrabItkImageMemory( tractsToDwiFilter->GetOutput() ); image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( parameters.m_SignalGen.GetGradientDirections() ) ); image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( parameters.m_SignalGen.m_Bvalue ) ); mitk::DiffusionPropertyHelper propertyHelper( image ); propertyHelper.InitializeImage(); if (file_extension=="") mitk::IOUtil::Save(image, "application/vnd.mitk.nii.gz", outName+".nii.gz"); else if (file_extension==".nii" || file_extension==".nii.gz") mitk::IOUtil::Save(image, "application/vnd.mitk.nii.gz", outName+file_extension); else mitk::IOUtil::Save(image, outName+file_extension); if (verbose) { std::vector< itk::TractsToDWIImageFilter< short >::ItkDoubleImgType::Pointer > volumeFractions = tractsToDwiFilter->GetVolumeFractions(); for (unsigned int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::IOUtil::Save(image, outName+"_Compartment"+boost::lexical_cast(k+1)+".nrrd"); } if (tractsToDwiFilter->GetPhaseImage().IsNotNull()) { mitk::Image::Pointer image = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkPhase = tractsToDwiFilter->GetPhaseImage(); image = mitk::GrabItkImageMemory( itkPhase.GetPointer() ); mitk::IOUtil::Save(image, outName+"_Phase.nrrd"); } if (tractsToDwiFilter->GetKspaceImage().IsNotNull()) { mitk::Image::Pointer image = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkImage = tractsToDwiFilter->GetKspaceImage(); image = mitk::GrabItkImageMemory( itkImage.GetPointer() ); mitk::IOUtil::Save(image, outName+"_kSpace.nrrd"); } } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt index 087a310ec5..22f786a84d 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt @@ -1,40 +1,42 @@ option(BUILD_DiffusionTractographyEvaluationCmdApps "Build commandline tools for diffusion fiber tractography evaluation" OFF) if(BUILD_DiffusionTractographyEvaluationCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionTractographyEvaluationcmdapps PeaksAngularError^^MitkFiberTracking TractometerMetrics^^MitkFiberTracking + TractPlausibility^^MitkFiberTracking + FitFibersToImage^^MitkFiberTracking # LocalDirectionalFiberPlausibility^^MitkFiberTracking # HAS TO USE NEW PEAK IMAGE FORMAT ) foreach(diffusionTractographyEvaluationcmdapp ${diffusionTractographyEvaluationcmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionTractographyEvaluationcmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} PACKAGE_DEPENDS ITK ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/FitFibersToImage.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/FitFibersToImage.cpp new file mode 100755 index 0000000000..c49dd60fb1 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/FitFibersToImage.cpp @@ -0,0 +1,450 @@ +/*=================================================================== + +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 + +using namespace std; +typedef itksys::SystemTools ist; +typedef itk::Image ItkUcharImgType; +typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType; +typedef mitk::DiffusionPropertyHelper DPH; +typedef itk::Point PointType; +typedef mitk::StickModel<> ModelType; +typedef mitk::BallModel<> BallModelType; + +std::vector SolveLinear(mitk::FiberBundle* inputTractogram, mitk::Image::Pointer inputImage, ModelType signalModel, BallModelType ) +{ + unsigned int num_gradients = signalModel.GetGradientList().size(); + DPH::ImageType::Pointer itkImage = DPH::GetItkVectorImage(inputImage); + + unsigned int* image_size = inputImage->GetDimensions(); + int sz_x = image_size[0]; + int sz_y = image_size[1]; + int sz_z = image_size[2]; + int num_voxels = sz_x*sz_y*sz_z; + + unsigned int num_unknowns = inputTractogram->GetNumFibers(); + unsigned int number_of_residuals = num_voxels * num_gradients; + + + // create linear system + MITK_INFO << "Num. unknowns: " << num_unknowns; + MITK_INFO << "Num. residuals: " << number_of_residuals; + + MITK_INFO << "Creating matrices ..."; + // Eigen::MatrixXf A = Eigen::MatrixXf::Constant(number_of_residuals, num_unknowns, 0); + // Eigen::VectorXf b = Eigen::VectorXf::Constant(number_of_residuals, 0); + + vigra::MultiArray<2, double> A(vigra::Shape2(number_of_residuals, num_unknowns), 0.0); + vigra::MultiArray<2, double> b(vigra::Shape2(number_of_residuals, 1), 0.0); + vigra::MultiArray<2, double> x(vigra::Shape2(num_unknowns, 1), 1.0); + + + itk::Image< double, 3>::Pointer temp_img1 = itk::Image< double, 3>::New(); + temp_img1->SetSpacing( itkImage->GetSpacing() ); + temp_img1->SetOrigin( itkImage->GetOrigin() ); + temp_img1->SetDirection( itkImage->GetDirection() ); + temp_img1->SetLargestPossibleRegion( itkImage->GetLargestPossibleRegion() ); + temp_img1->SetBufferedRegion( itkImage->GetBufferedRegion() ); + temp_img1->SetRequestedRegion( itkImage->GetRequestedRegion() ); + temp_img1->Allocate(); + temp_img1->FillBuffer(0.0); + + itk::Image< double, 3>::Pointer temp_img2 = itk::Image< double, 3>::New(); + temp_img2->SetSpacing( itkImage->GetSpacing() ); + temp_img2->SetOrigin( itkImage->GetOrigin() ); + temp_img2->SetDirection( itkImage->GetDirection() ); + temp_img2->SetLargestPossibleRegion( itkImage->GetLargestPossibleRegion() ); + temp_img2->SetBufferedRegion( itkImage->GetBufferedRegion() ); + temp_img2->SetRequestedRegion( itkImage->GetRequestedRegion() ); + temp_img2->Allocate(); + temp_img2->FillBuffer(0.0); + + MITK_INFO << "Filling matrices ..."; + unsigned int point_counter = 0; + vtkSmartPointer polydata = inputTractogram->GetFiberPolyData(); + for (int i=0; iGetNumFibers(); ++i) + { + vtkCell* cell = polydata->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + if (numPoints<2) + MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; + + for (int j=0; jGetPoint(j); + PointType p; + p[0]=p1[0]; + p[1]=p1[1]; + p[2]=p1[2]; + + DPH::ImageType::IndexType idx; + itkImage->TransformPhysicalPointToIndex(p, idx); + if (!itkImage->GetLargestPossibleRegion().IsInside(idx)) + continue; + + double* p2 = points->GetPoint(j+1); + ModelType::GradientType d; + d[0] = p[0]-p2[0]; + d[1] = p[1]-p2[1]; + d[2] = p[2]-p2[2]; + signalModel.SetFiberDirection(d); + ModelType::PixelType model_signal = signalModel.SimulateMeasurement(); + DPH::ImageType::PixelType measured_signal = itkImage->GetPixel(idx); + + int x = idx[0]; + int y = idx[1]; + int z = idx[2]; + + unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z; + + for (unsigned int k=0; kSetPixel(idx, temp_img1->GetPixel(idx) + (double)model_signal[k]); + A(num_gradients*linear_index + k, i) = temp_img1->GetPixel(idx); + +// temp_img2->SetPixel(idx, A(num_gradients*linear_index + k, i) ); + } + + point_counter++; + } + } + +// typedef itk::ImageFileWriter< itk::Image< double, 3> > WriterType; +// WriterType::Pointer writer = WriterType::New(); +// writer->SetFileName("/home/neher/Projects/TractPlausibility/model_signal.nrrd"); +// writer->SetInput(temp_img2); +// writer->Update(); + +// writer->SetFileName("/home/neher/Projects/TractPlausibility/measured_signal.nrrd"); +// writer->SetInput(temp_img1); +// writer->Update(); + + MITK_INFO << "Solving linear system"; + vigra::linalg::nonnegativeLeastSquares(A, b, x); + + std::vector weights; + for (unsigned int i=0; i SolveEvo(mitk::FiberBundle* inputTractogram, mitk::Image::Pointer inputImage, ModelType signalModel, BallModelType ballModel, int num_iterations=1000, double step=0.1) +{ + std::vector out_weights; + DPH::ImageType::Pointer itkImage = DPH::GetItkVectorImage(inputImage); + itk::VectorImage< double, 3>::Pointer simulatedImage = itk::VectorImage< double, 3>::New(); + simulatedImage->SetSpacing(itkImage->GetSpacing()); + simulatedImage->SetOrigin(itkImage->GetOrigin()); + simulatedImage->SetDirection(itkImage->GetDirection()); + simulatedImage->SetRegions(itkImage->GetLargestPossibleRegion()); + simulatedImage->SetVectorLength(itkImage->GetVectorLength()); + simulatedImage->Allocate(); + DPH::ImageType::PixelType zero_signal; + zero_signal.SetSize(itkImage->GetVectorLength()); + zero_signal.Fill(0); + simulatedImage->FillBuffer(zero_signal); + + vtkSmartPointer polydata = inputTractogram->GetFiberPolyData(); + + MITK_INFO << "INITIALIZING"; + std::vector< std::vector< ModelType::PixelType > > fiber_model_signals; + std::vector< std::vector< DPH::ImageType::IndexType > > fiber_image_indices; + + std::vector< int > fiber_indices; + int f = 0; + for (int i=0; iGetNumFibers(); ++i) + { + vtkCell* cell = polydata->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + +// if (numPoints<2) +// continue; + + std::vector< ModelType::PixelType > model_signals; + std::vector< DPH::ImageType::IndexType > image_indices; + + for (int j=0; jGetPoint(j); + PointType p; + p[0]=p1[0]; + p[1]=p1[1]; + p[2]=p1[2]; + + DPH::ImageType::IndexType idx; + itkImage->TransformPhysicalPointToIndex(p, idx); + if (!itkImage->GetLargestPossibleRegion().IsInside(idx)) + continue; + + double* p2 = points->GetPoint(j+1); + ModelType::GradientType d; + d[0] = p[0]-p2[0]; + d[1] = p[1]-p2[1]; + d[2] = p[2]-p2[2]; + signalModel.SetFiberDirection(d); + + model_signals.push_back(step*signalModel.SimulateMeasurement()); + image_indices.push_back(idx); + } + + fiber_model_signals.push_back(model_signals); + fiber_image_indices.push_back(image_indices); + fiber_indices.push_back(f); + out_weights.push_back(0); + ++f; + } + + BallModelType::PixelType ballSignal = step/10*ballModel.SimulateMeasurement(); + + double T_start = 0.0001; + double T_end = 0.000001; + double alpha = log(T_end/T_start); + + unsigned int num_gradients = signalModel.GetGradientList().size(); + for (int i=0; i model_signals = fiber_model_signals.at(f); + std::vector< DPH::ImageType::IndexType > image_indices = fiber_image_indices.at(f); + + double E_old = 0; + double E_new = 0; + + int add = std::rand()%2; +// add = 1; + int use_ball = std::rand()%2; + + if (add==0 && use_ball==0 && out_weights[f]GetPixel(idx); + ModelType::PixelType sVal = simulatedImage->GetPixel(idx); + ModelType::PixelType dVal = model_signals.at(c); + if (use_ball==1) + dVal = ballSignal; + + for (unsigned int g=0; g (rand()) / static_cast (RAND_MAX); + if (pGetPixel(idx); + ModelType::PixelType dVal = model_signals.at(c); + if (use_ball==1) + dVal = ballSignal; + + if (add==1) + sVal += dVal; + else + sVal -= dVal; + + simulatedImage->SetPixel(idx, sVal); + + ++c; + } + + if (add==1 && use_ball==0) + out_weights[f] += step; + else if (use_ball==0) + out_weights[f] -= step; + + } + + } + MITK_INFO << "Accepted: " << (float)accepted/fiber_indices.size(); + } + + typedef itk::ImageFileWriter< itk::VectorImage< double, 3> > WriterType; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName("/home/neher/Projects/TractPlausibility/model_signal.nrrd"); + writer->SetInput(simulatedImage); + writer->Update(); + + return out_weights; +} + + +/*! +\brief +*/ +int main(int argc, char* argv[]) +{ + mitkCommandLineParser parser; + + parser.setTitle("Fit Fibers To Image"); + parser.setCategory("Fiber Tracking Evaluation"); + parser.setDescription(""); + parser.setContributor("MIC"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input_tractogram", "i1", mitkCommandLineParser::InputFile, "Input Tractogram:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); + parser.addArgument("input_dMRI", "i2", mitkCommandLineParser::InputFile, "Input dMRI:", "input diffusion-weighted image", us::Any(), false); + parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); + parser.addArgument("step", "", mitkCommandLineParser::Float, "", "", us::Any()); + parser.addArgument("it", "", mitkCommandLineParser::Int, "", "", us::Any()); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + string fibFile = us::any_cast(parsedArgs["input_tractogram"]); + string dwiFile = us::any_cast(parsedArgs["input_dMRI"]); + string outRoot = us::any_cast(parsedArgs["out"]); + + int it = 1000; + if (parsedArgs.count("it")) + it = us::any_cast(parsedArgs["it"]); + + float step = 0.1; + if (parsedArgs.count("step")) + step = us::any_cast(parsedArgs["step"]); + + + + typedef DPH::GradientDirectionsContainerType GradientContainerType; + + try + { + mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); + + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "Fiberbundles"}, {}); + mitk::Image::Pointer inputImage = dynamic_cast(mitk::IOUtil::Load(dwiFile, &functor)[0].GetPointer()); + if (!DPH::IsDiffusionWeightedImage(inputImage)) + return EXIT_FAILURE; + + // resample fibers + float minSpacing = 1; + if(inputImage->GetGeometry()->GetSpacing()[0]GetGeometry()->GetSpacing()[1] && inputImage->GetGeometry()->GetSpacing()[0]GetGeometry()->GetSpacing()[2]) + minSpacing = inputImage->GetGeometry()->GetSpacing()[0]; + else if (inputImage->GetGeometry()->GetSpacing()[1] < inputImage->GetGeometry()->GetSpacing()[2]) + minSpacing = inputImage->GetGeometry()->GetSpacing()[1]; + else + minSpacing = inputImage->GetGeometry()->GetSpacing()[2]; + inputTractogram->ResampleLinear(minSpacing/10); + + // set up signal model + GradientContainerType* gradients = static_cast( inputImage->GetProperty(DPH::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); + + float b_value = 0; + inputImage->GetPropertyList()->GetFloatProperty(DPH::REFERENCEBVALUEPROPERTYNAME.c_str(), b_value); + + ModelType signalModel; + signalModel.SetDiffusivity(0.001); + signalModel.SetBvalue(b_value); + signalModel.SetGradientList(gradients); + + BallModelType ballModel; + ballModel.SetDiffusivity(0.001); + ballModel.SetBvalue(b_value); + ballModel.SetGradientList(gradients); + +// MITK_INFO << it << " " << step; +// std::vector weights = SolveLinear(inputTractogram, inputImage, signalModel, ballModel); + std::vector weights = SolveEvo(inputTractogram, inputImage, signalModel, ballModel, it, step); + + for (unsigned int i=0; iSetFiberWeight(i, weights.at(i)); + } + + mitk::IOUtil::Save(inputTractogram, outRoot + "fitted.fib"); + } + 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/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp new file mode 100755 index 0000000000..2b74fd5a94 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp @@ -0,0 +1,187 @@ +/*=================================================================== + +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 + +using namespace std; +typedef itksys::SystemTools ist; +typedef itk::Image ItkUcharImgType; +typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType; + +void CreateFolderStructure(const std::string& path) +{ + if (ist::PathExists(path)) + ist::RemoveADirectory(path); + + ist::MakeDirectory(path); + ist::MakeDirectory(path + "/known_tracts/"); + ist::MakeDirectory(path + "/plausible_tracts/"); + ist::MakeDirectory(path + "/implausible_tracts/"); +} + +ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) +{ + mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); + ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); + mitk::CastToItkImage(img, itkMask); + return itkMask; +} + +std::vector< MaskType > get_file_list(const std::string& path) +{ + std::vector< MaskType > mask_list; + + itk::Directory::Pointer dir = itk::Directory::New(); + + if (dir->Load(path.c_str())) + { + int n = dir->GetNumberOfFiles(); + for (int r = 0; r < n; r++) + { + const char *filename = dir->GetFile(r); + std::string ext = ist::GetFilenameExtension(filename); + if (ext==".nii" || ext==".nii.gz" || ext==".nrrd") + { + MaskType m(LoadItkMaskImage(path + '/' + filename), ist::GetFilenameWithoutExtension(filename)); + mask_list.push_back(m); + } + } + } + + return mask_list; +} + +mitk::FiberBundle::Pointer TransformToMRtrixSpace(mitk::FiberBundle::Pointer fib) +{ + mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); + vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); + matrix->Identity(); + matrix->SetElement(0,0,-matrix->GetElement(0,0)); + matrix->SetElement(1,1,-matrix->GetElement(1,1)); + geometry->SetIndexToWorldTransformByVtkMatrix(matrix); + + vtkSmartPointer transformFilter = vtkSmartPointer::New(); + transformFilter->SetInputData(fib->GetFiberPolyData()); + transformFilter->SetTransform(geometry->GetVtkTransform()); + transformFilter->Update(); + mitk::FiberBundle::Pointer transformed_fib = mitk::FiberBundle::New(transformFilter->GetOutput()); + return transformed_fib; +} + +/*! +\brief +*/ +int main(int argc, char* argv[]) +{ + mitkCommandLineParser parser; + + parser.setTitle("Tract Plausibility"); + parser.setCategory("Fiber Tracking Evaluation"); + parser.setDescription(""); + parser.setContributor("MIC"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); + parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); + parser.addArgument("reference_mask_folder", "m", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", false); + parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", false); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + string fibFile = us::any_cast(parsedArgs["input"]); + string gray_matter_mask = us::any_cast(parsedArgs["gray_matter_mask"]); + string reference_mask_folder = us::any_cast(parsedArgs["reference_mask_folder"]); + string outRoot = us::any_cast(parsedArgs["out"]); + + try + { + CreateFolderStructure(outRoot); + + std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder); + if (known_tract_masks.empty()) + return EXIT_FAILURE; + + ItkUcharImgType::Pointer gm_image = LoadItkMaskImage(gray_matter_mask); + + mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); + + // filter gray matter fibers + mitk::FiberBundle::Pointer not_gm_fibers = inputTractogram->ExtractFiberSubset(gm_image, false, true, true); + mitk::IOUtil::Save(not_gm_fibers, outRoot + "/implausible_tracts/no_gm_endings.trk"); + inputTractogram = inputTractogram->ExtractFiberSubset(gm_image, false, false, true); + + // resample fibers + float minSpacing = 1; + if(std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[1] && std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[2]) + minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]; + else if (std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] < std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]) + minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[1]; + else + minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]; + inputTractogram->ResampleLinear(minSpacing/5); + + // find known tracts via overlap match + mitk::FiberBundle::Pointer all_known_tracts = nullptr; + for ( MaskType mask : known_tract_masks ) + { + ItkUcharImgType::Pointer mask_image = std::get<0>(mask); + std::string mask_name = std::get<1>(mask); + mitk::FiberBundle::Pointer known_tract = inputTractogram->ExtractFiberSubset(mask_image, true, false, false, 0.9, false); + mitk::IOUtil::Save(known_tract, outRoot + "/known_tracts/" + mask_name + ".trk"); + + if (all_known_tracts.IsNull()) + all_known_tracts = mitk::FiberBundle::New(known_tract->GetFiberPolyData()); + else + all_known_tracts = all_known_tracts->AddBundle(known_tract); + } + mitk::IOUtil::Save(all_known_tracts, outRoot + "/known_tracts/all_known_tracts.trk"); + mitk::IOUtil::Save(TransformToMRtrixSpace(all_known_tracts), outRoot + "/known_tracts/all_known_tracts.fib"); + + mitk::FiberBundle::Pointer remaining_tracts = inputTractogram->SubtractBundle(all_known_tracts); + mitk::IOUtil::Save(remaining_tracts, outRoot + "/plausible_tracts/remaining_tracts.trk"); + mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), outRoot + "/plausible_tracts/remaining_tracts.fib"); + MITK_INFO << "step_size: " << minSpacing/5; + } + 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/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp index 5c897e61b1..3e44654736 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp @@ -1,1540 +1,1539 @@ /*=================================================================== 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 "QmitkFiberProcessingView.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include "mitkNodePredicateDataType.h" #include #include #include #include // ITK #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() : QmitkAbstractView() , m_Controls( 0 ) , m_CircleCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(1) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { RemoveObservers(); } void QmitkFiberProcessingView::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::QmitkFiberProcessingViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); m_Controls->label_17->setVisible(false); m_Controls->m_FiberExtractionFractionBox->setVisible(false); } UpdateGui(); } void QmitkFiberProcessingView::OnMaskExtractionChanged() { if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2) { m_Controls->label_17->setVisible(true); m_Controls->m_FiberExtractionFractionBox->setVisible(true); m_Controls->m_BothEnds->setVisible(false); } else { m_Controls->label_17->setVisible(false); m_Controls->m_FiberExtractionFractionBox->setVisible(false); if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) m_Controls->m_BothEnds->setVisible(true); } } void QmitkFiberProcessingView::SetFocus() { m_Controls->toolBoxx->setFocus(); } void QmitkFiberProcessingView::Modify() { switch (m_Controls->m_ModificationMethodBox->currentIndex()) { case 0: { ResampleSelectedBundlesSpline(); break; } case 1: { ResampleSelectedBundlesLinear(); break; } case 2: { CompressSelectedBundles(); break; } case 3: { DoImageColorCoding(); break; } case 4: { MirrorFibers(); break; } case 5: { WeightFibers(); break; } case 6: { DoCurvatureColorCoding(); break; } case 7: { DoWeightColorCoding(); break; } } } void QmitkFiberProcessingView::WeightFibers() { float weight = this->m_Controls->m_BundleWeightBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->SetFiberWeights(weight); } } void QmitkFiberProcessingView::Remove() { switch (m_Controls->m_RemovalMethodBox->currentIndex()) { case 0: { RemoveDir(); break; } case 1: { PruneBundle(); break; } case 2: { ApplyCurvatureThreshold(); break; } case 3: { RemoveWithMask(false); break; } case 4: { RemoveWithMask(true); break; } } } void QmitkFiberProcessingView::Extract() { switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { case 0: { ExtractWithPlanarFigure(); break; } case 1: { switch (m_Controls->m_ExtractionBoxMask->currentIndex()) { { case 0: ExtractWithMask(true, false); break; } { case 1: ExtractWithMask(true, true); break; } { case 2: ExtractWithMask(false, false); break; } { case 3: ExtractWithMask(false, true); break; } } break; } } } void QmitkFiberProcessingView::PruneBundle() { int minLength = this->m_Controls->m_PruneFibersMinBox->value(); int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (!fib->RemoveShortFibers(minLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); else if (!fib->RemoveLongFibers(maxLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyCurvatureThreshold() { int angle = this->m_Controls->m_CurvSpinBox->value(); int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(angle); filter->SetDistance(dist); filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); filter->Update(); mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByOrientation(); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveDir() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); vnl_vector_fixed dir; dir[0] = m_Controls->m_ExtractDirX->value(); dir[1] = m_Controls->m_ExtractDirY->value(); dir[2] = m_Controls->m_ExtractDirZ->value(); fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*M_PI/180)); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { if (m_MaskImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itkUCharImageType::Pointer mask = itkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } node->SetData(newFib); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert) { if (m_MaskImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_MaskImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); itkUCharImageType::Pointer mask = itkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->ExtractFiberSubset(mask, !onlyEnds, invert, m_Controls->m_BothEnds->isChecked(), m_Controls->m_FiberExtractionFractionBox->value()); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } DataNode::Pointer newNode = DataNode::New(); newNode->SetData(newFib); if (invert) { name += "_not"; if (onlyEnds) name += "-ending-in-mask"; else name += "-passing-mask"; } else { if (onlyEnds) name += "_ending-in-mask"; else name += "_passing-mask"; } newNode->SetName(name.toStdString()); GetDataStorage()->Add(newNode); node->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateRoiImage() { if (m_SelectedPF.empty()) return; mitk::BaseGeometry::Pointer geometry; if (!m_SelectedFB.empty()) { mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); geometry = fib->GetGeometry(); } else if (m_SelectedImage) geometry = m_SelectedImage->GetGeometry(); else return; itk::Vector spacing = geometry->GetSpacing(); spacing /= m_UpsamplingFactor; mitk::Point3D newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); itk::Matrix direction; itk::ImageRegion<3> imageRegion; for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); m_PlanarFigureImage = itkUCharImageType::New(); m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin m_PlanarFigureImage->SetDirection( direction ); // Set the image direction m_PlanarFigureImage->SetRegions( imageRegion ); m_PlanarFigureImage->Allocate(); m_PlanarFigureImage->FillBuffer( 0 ); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); std::string name = m_SelectedPF.at(0)->GetName(); WritePfToImage(m_SelectedPF.at(0), tmpImage); for (unsigned int i=1; iGetName(); WritePfToImage(m_SelectedPF.at(i), tmpImage); } DataNode::Pointer node = DataNode::New(); tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName(name); this->GetDataStorage()->Add(node); } void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node->GetData())) { m_PlanarFigure = dynamic_cast(node->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } else if (dynamic_cast(node->GetData())) { DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; iSize(); i++) { WritePfToImage(children->at(i), image); } } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = 0.5; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetExtentInMM(0) / spacing[0]; size[1] = planegeo->GetExtentInMM(1) / spacing[1]; size[2] = 1; resampler->SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; corrorig[2] += 0; planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(unsigned int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) sigma[d] = gausssigma * image->GetSpacing()[d]; double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, itkUCharImageType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". itkUCharImageType::Pointer newMaskImage = itkUCharImageType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 1 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // Create VTK polydata object of polyline contour vtkPoints *points = vtkPoints::New(); PlanarFigure::PolyLineType::const_iterator it; unsigned int numberOfPoints = 0; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected image Point2D point2D = *it; planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); planarFigurePlaneGeometry->Map( point2D, point3D ); // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { float bounds[2] = {0,0}; bounds[0] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); bounds[1] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); imageGeometry3D->WorldToIndex( point3D, point3D ); if (point3D[i0]<0) point3D[i0] = 0.0; else if (point3D[i0]>bounds[0]) point3D[i0] = bounds[0]-0.001; if (point3D[i1]<0) point3D[i1] = 0.0; else if (point3D[i1]>bounds[1]) point3D[i1] = bounds[1]-0.001; points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } else { imageGeometry3D->WorldToIndex( point3D, point3D ); // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } } polyline->SetPoints( points ); points->Delete(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) ptIds[i] = i; polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInputData( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< itkUCharImageType > ImageImportType; typedef itk::VTKImageExport< itkUCharImageType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; typename ImageType::SizeType uppersize = {{0,0,0}}; while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) itimage.Set(0); else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); // crop internal mask typedef itk::RegionOfInterestImageFilter< itkUCharImageType, itkUCharImageType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(cropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); Image::Pointer tmpImage2 = Image::New(); tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 ); typedef itk::ImageRegionIteratorWithIndex IteratorType; IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); imageIterator.GoToBegin(); while ( !imageIterator.IsAtEnd() ) { unsigned char val = imageIterator.Value(); if (val>0) { itk::Index<3> index = imageIterator.GetIndex(); Point3D point; point[0] = index[0]; point[1] = index[1]; point[2] = index[2]; intImageGeometry3D->IndexToWorld(point, point); pfImageGeometry3D->WorldToIndex(point, point); point[i0] += 0.5; point[i1] += 0.5; index[0] = point[0]; index[1] = point[1]; index[2] = point[2]; if (pfImageGeometry3D->IsIndexInside(index)) m_PlanarFigureImage->SetPixel(index, 1); } ++imageIterator; } // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void QmitkFiberProcessingView::UpdateGui() { m_Controls->m_FibLabel->setText("mandatory"); m_Controls->m_PfLabel->setText("needed for extraction"); m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_RemoveButton->setEnabled(false); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); m_Controls->m_ExtractFibersButton->setEnabled(false); m_Controls->m_ModifyButton->setEnabled(false); m_Controls->m_CopyBundle->setEnabled(false); m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); // disable alle frames m_Controls->m_BundleWeightFrame->setVisible(false); m_Controls->m_ExtactionFramePF->setVisible(false); m_Controls->m_RemoveDirectionFrame->setVisible(false); m_Controls->m_RemoveLengthFrame->setVisible(false); m_Controls->m_RemoveCurvatureFrame->setVisible(false); m_Controls->m_SmoothFibersFrame->setVisible(false); m_Controls->m_CompressFibersFrame->setVisible(false); m_Controls->m_ColorFibersFrame->setVisible(false); m_Controls->m_MirrorFibersFrame->setVisible(false); m_Controls->m_MaskExtractionFrame->setVisible(false); m_Controls->m_ColorMapBox->setVisible(false); bool pfSelected = !m_SelectedPF.empty(); bool fibSelected = !m_SelectedFB.empty(); bool multipleFibsSelected = (m_SelectedFB.size()>1); bool maskSelected = m_MaskImageNode.IsNotNull(); bool imageSelected = m_SelectedImage.IsNotNull(); // toggle visibility of elements according to selected method switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) { case 0: m_Controls->m_ExtactionFramePF->setVisible(true); break; case 1: m_Controls->m_MaskExtractionFrame->setVisible(true); break; } switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) { case 0: m_Controls->m_RemoveDirectionFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 1: m_Controls->m_RemoveLengthFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 2: m_Controls->m_RemoveCurvatureFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 3: break; case 4: break; } switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) { case 0: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 1: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 2: m_Controls->m_CompressFibersFrame->setVisible(true); break; case 3: m_Controls->m_ColorFibersFrame->setVisible(true); m_Controls->m_ColorMapBox->setVisible(true); break; case 4: m_Controls->m_MirrorFibersFrame->setVisible(true); if (m_SelectedSurfaces.size()>0) m_Controls->m_ModifyButton->setEnabled(true); break; case 5: m_Controls->m_BundleWeightFrame->setVisible(true); break; case 6: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 7: m_Controls->m_ColorFibersFrame->setVisible(true); break; } // are fiber bundles selected? if ( fibSelected ) { m_Controls->m_CopyBundle->setEnabled(true); m_Controls->m_ModifyButton->setEnabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); // one bundle and one planar figure needed to extract fibers if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } // more than two bundles needed to join/subtract if (multipleFibsSelected) { m_Controls->m_FibLabel->setText("multiple bundles selected"); m_Controls->m_JoinBundles->setEnabled(true); m_Controls->m_SubstractBundles->setEnabled(true); } if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str())); m_Controls->m_RemoveButton->setEnabled(true); } } // are planar figures selected? if (pfSelected) { if ( fibSelected || m_SelectedImage.IsNotNull()) m_Controls->m_GenerateRoiImage->setEnabled(true); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); } else m_Controls->PFCompoNOTButton->setEnabled(true); } // is image selected if (imageSelected || maskSelected) { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); } } void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node ) { for (auto fnode: m_SelectedFB) if (node == fnode) { m_SelectedFB.clear(); break; } berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* ) { if (!m_Controls->m_InteractiveBox->isChecked()) { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } } void QmitkFiberProcessingView::OnEndInteraction() { if (m_Controls->m_InteractiveBox->isChecked()) ExtractWithPlanarFigure(true); } void QmitkFiberProcessingView::AddObservers() { typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType; for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) { figure->RemoveAllObservers(); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction); m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); } } } void QmitkFiberProcessingView::RemoveObservers() { for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) figure->RemoveAllObservers(); } } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { RemoveObservers(); //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection std::vector lastSelectedFB = m_SelectedFB; m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; m_MaskImageNode = nullptr; for (auto node: nodes) { if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back(node); else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) { m_SelectedImage = dynamic_cast(node->GetData()); bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) m_MaskImageNode = node; } else if (dynamic_cast(node->GetData())) m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked()) m_SelectedFB = lastSelectedFB; // if no fibers or surfaces are selected, select topmost if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedFB.clear(); m_SelectedFB.push_back(nodes->at(i)); } } } // if no plar figure is selected, select topmost if (m_SelectedPF.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedPF.clear(); m_SelectedPF.push_back(nodes->at(i)); } } } AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); } void QmitkFiberProcessingView::OnDrawCircle() { mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetPlaneGeometry( emptygeometry ); //set desired data to DataNode where Planarfigure is stored mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); newNode->SetBoolProperty("planarfigure.3drendering", true); newNode->SetBoolProperty("planarfigure.3drendering.fill", true); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode(newNode); } // figure drawn on the topmost layer / image GetDataStorage()->Add(newNode ); RemoveObservers(); for(unsigned int i = 0; i < m_SelectedPF.size(); i++) m_SelectedPF[i]->SetSelected(false); newNode->SetSelected(true); m_SelectedPF.clear(); m_SelectedPF.push_back(newNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive) { if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); return; } std::vector fiberBundles = m_SelectedFB; mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); if (interactive && m_Controls->m_InteractiveBox->isChecked()) { - if (extFB->GetNumFibers()<=0) - continue; - if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } float op = 5.0/sqrt(fib->GetNumFibers()); + float currentOp = 0; + fiberBundles.at(i)->GetFloatProperty("opacity", currentOp); - fib->SetFiberColors(255, 255, 255); - fiberBundles.at(i)->SetFloatProperty("opacity", op); - fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); + if (currentOp!=op) + { + fib->SetFiberColors(255, 255, 255); + fiberBundles.at(i)->SetFloatProperty("opacity", op); + fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); + } m_InteractiveNode->SetData(extFB); - - if (auto renderWindowPart = this->GetRenderWindowPart()) - renderWindowPart->RequestUpdate(); } else { if (extFB->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); QString name(fiberBundles.at(i)->GetName().c_str()); name += "*"; node->SetName(name.toStdString()); fiberBundles.at(i)->SetVisibility(false); GetDataStorage()->Add(node); } } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("AND"); newPFCNode->SetData(PFCAnd); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); RemoveObservers(); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("OR"); newPFCNode->SetData(PFCOr); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("NOT"); newPFCNode->SetData(PFCNot); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode ) { pfc->SetSelected(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(pfc, parentNode); else GetDataStorage()->Add(pfc); for (auto child : children) { if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); newChild->SetBoolProperty("planarfigure.3drendering", true); newChild->SetBoolProperty("planarfigure.3drendering.fill", true); GetDataStorage()->Add(newChild, pfc); GetDataStorage()->Remove(child); } else if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); std::vector< mitk::DataNode::Pointer > grandChildVector; mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) grandChildVector.push_back(*it); AddCompositeToDatastorage(newChild, grandChildVector, pfc); GetDataStorage()->Remove(child); } } UpdateGui(); } void QmitkFiberProcessingView::CopyBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); node->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newFib); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } UpdateGui(); } void QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iAddBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); name += "+"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::SubstractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); if (newBundle.IsNull()) break; name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } if (newBundle.IsNull()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); return; } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::ResampleSelectedBundlesSpline() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleSpline(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ResampleSelectedBundlesLinear() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleLinear(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::CompressSelectedBundles() { double factor = this->m_Controls->m_ErrorThresholdBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->Compress(factor); fib->ColorFibersByOrientation(); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoCurvatureColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoWeightColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::MirrorFibers() { unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (m_SelectedImage.IsNotNull()) fib->SetReferenceGeometry(m_SelectedImage->GetGeometry()); fib->MirrorFibers(axis); } for (auto surf : m_SelectedSurfaces) { vtkSmartPointer poly = surf->GetVtkPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); for (int i=0; iGetNumberOfPoints(); i++) { double* point = poly->GetPoint(i); point[axis] *= -1; vtkNewPoints->InsertNextPoint(point); } poly->SetPoints(vtkNewPoints); surf->CalculateBoundingBox(); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } }