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 <mitkITKImageImport.h>
#include <mitkImageCast.h>
#include <itkImageRegionIterator.h>
#include <mitkProperties.h>
#include <vnl/algo/vnl_matrix_inverse.h>
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<mitk::GradientDirectionsProperty*>( 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<std::vector<int> > dirIndices;
for(GradientDirectionsContainerType::ConstIterator gdcitNew = newDirs->Begin();
gdcitNew != newDirs->End(); ++gdcitNew)
{
dirIndices.push_back(std::vector<int>(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; i<dirIndices.size(); i++)
{
// do the averaging
const unsigned int numavg = dirIndices[i].size();
unsigned int sum = 0;
for(unsigned int j=0; j<numavg; j++)
{
//MITK_INFO << newVec[i] << " << " << oldVec[dirIndices[i].at(j)];
sum += oldVec[dirIndices[i].at(j)];
}
if(numavg == 0)
{
MITK_ERROR << "VectorImage: Error on averaging. Possibly due to corrupted data";
return;
}
newVec[i] = sum / numavg;
}
newIt.Set(newVec);
++newIt;
++oldIt;
}
mitk::GrabItkImageMemory( newITKImage, m_Image );
m_Image->GetPropertyList()->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<mitk::GradientDirectionsProperty*>( 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<double> 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<mitk::GradientDirectionsProperty*>( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer())->GetGradientDirectionsContainer();
MeasurementFrameType measurementFrame = GetMeasurementFrame(m_Image);
measurementFrame = vnl_matrix_inverse<double>(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<double>(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<double> 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<mitk::GradientDirectionsProperty*>( 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<double> 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<mitk::BValueMapProperty*>(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<mitk::GradientDirectionsProperty*>( 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<mitk::GradientDirectionsProperty*>( m_Image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer();
float b_value = static_cast<mitk::FloatProperty*>(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<mitk::MeasurementFrameProperty *>( 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<mitk::Image *>(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<mitk::FloatProperty *>(image->GetProperty(REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer());
if( referenceBValue.IsNull() )
{
isDiffusionWeightedImage = false;
}
}
unsigned int gradientDirections( 0 );
if( isDiffusionWeightedImage )
{
mitk::GradientDirectionsProperty::Pointer gradientDirectionsProperty = dynamic_cast<mitk::GradientDirectionsProperty *>(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<mitk::BValueMapProperty *>(image->GetProperty(BVALUEMAPPROPERTYNAME.c_str()).GetPointer())->GetBValueMap();
}
float mitk::DiffusionPropertyHelper::GetReferenceBValue(const mitk::Image *image)
{
return dynamic_cast<mitk::FloatProperty *>(image->GetProperty(REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer())->GetValue();
}
const mitk::DiffusionPropertyHelper::MeasurementFrameType & mitk::DiffusionPropertyHelper::GetMeasurementFrame(const mitk::Image *image)
{
mitk::MeasurementFrameProperty::Pointer mf = dynamic_cast<mitk::MeasurementFrameProperty *>( 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<mitk::GradientDirectionsProperty *>(image->GetProperty(ORIGINALGRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer())->GetGradientDirectionsContainer();
}
mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer mitk::DiffusionPropertyHelper::GetGradientContainer(const mitk::Image *image)
{
return dynamic_cast<mitk::GradientDirectionsProperty *>(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 <itkMetaDataObject.h>
#include <vtkPolyData.h>
#include <vtkDataReader.h>
#include <vtkPolyDataReader.h>
#include <vtkMatrix4x4.h>
#include <vtkPolyLine.h>
#include <vtkCellArray.h>
#include <vtkDataArray.h>
#include <vtkFloatArray.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <itksys/SystemTools.hxx>
#include <tinyxml.h>
#include <vtkCleanPolyData.h>
#include <mitkTrackvis.h>
#include <mitkCustomMimeType.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkXMLDataReader.h>
#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<itk::SmartPointer<mitk::BaseData> > mitk::FiberBundleVtkReader::Read()
{
std::vector<itk::SmartPointer<mitk::BaseData> > 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<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
reader->SetFileName( this->GetInputLocation().c_str() );
if (reader->IsFilePolyData())
{
reader->Update();
if ( reader->GetOutput() != nullptr )
{
vtkSmartPointer<vtkPolyData> fiberPolyData = reader->GetOutput();
FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData);
vtkSmartPointer<vtkFloatArray> weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS"));
if (weights!=nullptr)
{
// float weight=0;
- for (int i=0; i<weights->GetSize(); i++)
+ for (int i=0; i<weights->GetNumberOfValues(); 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<vtkUnsignedCharArray> 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<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
reader->SetFileName( this->GetInputLocation().c_str() );
if ( reader->CanReadFile(this->GetInputLocation().c_str()) )
{
reader->Update();
if ( reader->GetOutput() != nullptr )
{
vtkSmartPointer<vtkPolyData> fiberPolyData = reader->GetOutput();
FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData);
vtkSmartPointer<vtkFloatArray> weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS"));
if (weights!=nullptr)
{
// float weight=0;
// for (int i=0; i<weights->GetSize(); i++)
// if (!mitk::Equal(weights->GetValue(i),weight,0.00001))
// {
// MITK_INFO << "Weight: " << weights->GetValue(i);
// weight = weights->GetValue(i);
// }
fiberBundle->SetFiberWeights(weights);
}
vtkSmartPointer<vtkUnsignedCharArray> 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<vtkPolyData> fiberPolyData = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkCellArray> cellArray = vtkSmartPointer<vtkCellArray>::New();
vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::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<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::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<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::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 <mitkDiffusionSignalModel.h>
#include <itkPointShell.h>
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<OtherType>* 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<double> sticks;
switch (order)
{
case 1:
m_NumSticks = 12;
sticks = itk::PointShell<12, vnl_matrix_fixed<double, 3, 12> >::DistributePointShell()->as_matrix();
break;
case 2:
m_NumSticks = 42;
sticks = itk::PointShell<42, vnl_matrix_fixed<double, 3, 42> >::DistributePointShell()->as_matrix();
break;
case 3:
m_NumSticks = 92;
sticks = itk::PointShell<92, vnl_matrix_fixed<double, 3, 92> >::DistributePointShell()->as_matrix();
break;
case 4:
m_NumSticks = 162;
sticks = itk::PointShell<162, vnl_matrix_fixed<double, 3, 162> >::DistributePointShell()->as_matrix();
break;
case 5:
m_NumSticks = 252;
sticks = itk::PointShell<252, vnl_matrix_fixed<double, 3, 252> >::DistributePointShell()->as_matrix();
break;
default:
m_NumSticks = 42;
sticks = itk::PointShell<42, vnl_matrix_fixed<double, 3, 42> >::DistributePointShell()->as_matrix();
break;
}
for (int i=0; i<m_NumSticks; i++)
{
GradientType stick;
stick[0] = sticks.get(0,i); stick[1] = sticks.get(1,i); stick[2] = sticks.get(2,i);
stick.Normalize();
m_Sticks.push_back(stick);
}
}
unsigned int GetNumSticks(){ return m_NumSticks; }
GradientListType GetSticks(){ return m_Sticks; }
protected:
GradientType GetRandomDirection();
double m_BValue; ///< b-value used to generate the artificial signal
double m_Diffusivity; ///< Scalar diffusion constant
GradientListType m_Sticks; ///< Stick container
unsigned int m_NumSticks; ///< Number of sticks
bool m_RandomizeSticks;
};
}
#include "mitkAstroStickModel.cpp"
#endif
diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h
index 98a0bd7e33..3629525dd3 100644
--- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h
+++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/SignalModels/mitkBallModel.h
@@ -1,77 +1,76 @@
/*===================================================================
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_BallModel_H
#define _MITK_BallModel_H
#include <mitkDiffusionSignalModel.h>
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<OtherType>* 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 <MitkFiberTrackingExports.h>
#include <itkVariableLengthVector.h>
#include <itkVector.h>
#include <itkImage.h>
#include <itkMersenneTwisterRandomVariateGenerator.h>
#include <vnl/vnl_vector_fixed.h>
+#include <mitkDiffusionPropertyHelper.h>
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<double, 3> ItkDoubleImgType;
typedef itk::VariableLengthVector< ScalarType > PixelType;
typedef itk::Vector<double,3> GradientType;
typedef std::vector<GradientType> 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; i<gradients->Size(); ++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 <mitkDiffusionSignalModel.h>
namespace mitk {
/**
* \brief Generates constant direction independent signal.
*
*/
template< class ScalarType = double >
class DotModel : public DiffusionSignalModel< ScalarType >
{
public:
DotModel();
template< class OtherType >DotModel(DotModel<OtherType>* 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 <mitkDiffusionSignalModel.h>
#include <mitkImage.h>
#include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
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<OtherType>* 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<short,short,float,2,QBALL_ODFSIZE> 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<double> 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<double> 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 <mitkDiffusionSignalModel.h>
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<OtherType>* 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 <mitkDiffusionSignalModel.h>
#include <itkDiffusionTensor3D.h>
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<OtherType>* 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<double, 3, 3> 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<double, 3, 3> 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 <mitkPlanarCircle.h>
#include <mitkPlanarPolygon.h>
#include <mitkPlanarFigureComposite.h>
#include "mitkImagePixelReadAccessor.h"
#include <mitkPixelTypeMultiplex.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkUnsignedCharArray.h>
#include <vtkPolyLine.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkIdFilter.h>
#include <vtkClipPolyData.h>
#include <vtkPlane.h>
#include <vtkDoubleArray.h>
#include <vtkKochanekSpline.h>
#include <vtkParametricFunctionSource.h>
#include <vtkParametricSpline.h>
#include <vtkPolygon.h>
#include <vtkCleanPolyData.h>
#include <cmath>
#include <boost/progress.hpp>
#include <vtkTransformPolyDataFilter.h>
#include <mitkTransferFunction.h>
#include <vtkLookupTable.h>
#include <mitkLookupTable.h>
#include <vtkCardinalSpline.h>
const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs";
using namespace std;
mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData )
: m_NumFibers(0)
{
m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
m_FiberWeights->SetName("FIBER_WEIGHTS");
m_FiberPolyData = vtkSmartPointer<vtkPolyData>::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<vtkPolyData> mitk::FiberBundle::GeneratePolyDataByIds(std::vector<long> fiberIds)
{
vtkSmartPointer<vtkPolyData> newFiberPolyData = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkCellArray> newLineSet = vtkSmartPointer<vtkCellArray>::New();
vtkSmartPointer<vtkPoints> newPointSet = vtkSmartPointer<vtkPoints>::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<vtkCell> fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber);
vtkSmartPointer<vtkPoints> fibPoints = fiber->GetPoints();
vtkSmartPointer<vtkPolyLine> newFiber = vtkSmartPointer<vtkPolyLine>::New();
newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() );
for(int i=0; i<fibPoints->GetNumberOfPoints(); 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<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
// add current fiber bundle
vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers());
unsigned int counter = 0;
for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double p[3];
points->GetPoint(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; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
{
vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double p[3];
points->GetPoint(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<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
std::vector< std::vector< itk::Point<float, 3> > > points1;
for( int i=0; i<m_NumFibers; i++ )
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
if (points==nullptr || numPoints<=0)
continue;
itk::Point<float, 3> start = GetItkPoint(points->GetPoint(0));
itk::Point<float, 3> end = GetItkPoint(points->GetPoint(numPoints-1));
points1.push_back( {start, end} );
}
std::vector< std::vector< itk::Point<float, 3> > > points2;
for( int i=0; i<fib->GetNumFibers(); i++ )
{
vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
if (points==nullptr || numPoints<=0)
continue;
itk::Point<float, 3> start = GetItkPoint(points->GetPoint(0));
itk::Point<float, 3> 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; j<points2.size(); j++)
{
auto v1 = points1.at(i);
auto v2 = points2.at(j);
unsigned int matches = 0;
unsigned int reverse_matches = 0;
for (unsigned int c=0; c<v1.size(); c++)
{
if (v1[c].SquaredEuclideanDistanceTo(v2[c])<mitk::eps)
matches++;
if (v1[v1.size() - c - 1].SquaredEuclideanDistanceTo(v2[c])<mitk::eps)
reverse_matches++;
}
if (matches==v1.size() || reverse_matches==v1.size())
{
match = true;
j=points2.size();
}
}
#pragma omp critical
if (!match)
ids.push_back(i);
}
for( int i : ids )
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
if (points==nullptr || numPoints<=0)
continue;
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for( int j=0; j<numPoints; j++)
{
vtkIdType id = vNewPoints->InsertNextPoint(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<float, 3> mitk::FiberBundle::GetItkPoint(double point[3])
{
itk::Point<float, 3> 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<vtkPolyData> fiberPD, bool updateGeometry)
{
if (fiberPD == nullptr)
this->m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
else
m_FiberPolyData->DeepCopy(fiberPD);
m_NumFibers = m_FiberPolyData->GetNumberOfLines();
if (updateGeometry)
UpdateFiberGeometry();
GenerateFiberIds();
ColorFibersByOrientation();
}
/*
* return vtkPolyData
*/
vtkSmartPointer<vtkPolyData> 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<vtkUnsignedCharArray>::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; fi<numOfFibers; ++fi) {
vtkIdType* idList; // contains the point id's of the line
vtkIdType pointsPerFiber; // number of points for current line
fiberList->GetNextCell(pointsPerFiber, idList);
/* single fiber checkpoints: is number of points valid */
if (pointsPerFiber > 1)
{
/* operate on points of single fiber */
for (int i=0; i <pointsPerFiber; ++i)
{
/* process all points elastV[0]ept starting and endpoint for calculating color value take current point, previous point and next point */
if (i<pointsPerFiber-1 && 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<vtkUnsignedCharArray>::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<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::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; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
++disp;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
// calculate curvatures
for (int j=0; j<numPoints; j++)
{
double dist = 0;
int c = j;
std::vector< vnl_vector_fixed< float, 3 > > vectors;
vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0);
while(dist<window/2 && c>1)
{
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(dist<window/2 && c<numPoints-1)
{
double p1[3];
points->GetPoint(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; c<vectors.size(); c++)
{
double angle = dot_product(meanV, vectors.at(c));
if (angle>1.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 (dev<min)
min = dev;
if (dev>max)
max = dev;
}
}
unsigned int count = 0;
for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
for (int j=0; j<numPoints; j++)
{
double color[3];
double dev = values.at(count);
if (normalize)
dev = (dev-min)/(max-min);
else if (dev>1)
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; i<m_FiberColors->GetNumberOfTuples(); 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; i<m_FiberColors->GetNumberOfTuples(); 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 <typename TPixel>
void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize)
{
m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
m_FiberColors->SetNumberOfComponents(4);
m_FiberColors->SetName("FIBER_COLORS");
mitk::ImagePixelReadAccessor<TPixel,3> 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<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::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; i<m_FiberPolyData->GetNumberOfPoints(); ++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 (pixelValue<min)
min = pixelValue;
}
for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++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<vtkUnsignedCharArray>::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<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::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; i<m_NumFibers; i++)
{
double weight = this->GetFiberWeight(i);
if (weight>max)
max = weight;
if (weight<min)
min = weight;
}
if (fabs(max-min)<0.00001)
{
max = 1;
min = 0;
}
for (int i=0; i<m_NumFibers; i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
double weight = this->GetFiberWeight(i);
for (int j=0; j<numPoints; j++)
{
float v = weight;
if (normalize)
v = (v-min)/(max-min);
else if (v>1)
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<vtkUnsignedCharArray>::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; i<m_FiberPolyData->GetNumberOfPoints(); ++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<vtkIdFilter> idFiberFilter = vtkSmartPointer<vtkIdFilter>::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<vtkPolyData> PolyData = m_FiberPolyData;
- if (anyPoint)
+ if (anyPoint && do_resampling)
{
float minSpacing = 1;
if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
MITK_INFO << "Extracting fibers with mask image";
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp;
vtkCell* cell = PolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkCell* cellOriginal = m_FiberPolyData->GetCell(i);
int numPointsOriginal = cellOriginal->GetNumberOfPoints();
vtkPoints* pointsOriginal = cellOriginal->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
if (numPoints>1 && numPointsOriginal)
{
if (anyPoint)
{
int inside = 0;
int outside = 0;
if (!invert)
{
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
itk::Point<float, 3> 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; k<numPoints; k++)
{
double* p = points->GetPoint(k);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
else
{
bool includeFiber = true;
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
itk::Point<float, 3> 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; k<numPoints; k++)
{
double* p = points->GetPoint(k);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
}
else
{
double* start = pointsOriginal->GetPoint(0);
itk::Point<float, 3> 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<float, 3> 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; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 )
{
for (int j=0; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(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; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(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; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
}
}
vtkNewCells->InsertNextCell(container);
}
if (vtkNewCells->GetNumberOfCells()<=0)
return nullptr;
vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::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]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->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<vtkPolyData> PolyData =fibCopy->GetFiberPolyData();
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
MITK_INFO << "Cutting fibers";
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp;
vtkCell* cell = PolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
if (numPoints>1)
{
int newNumPoints = 0;
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
itk::Point<float, 3> 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<vtkPolyLine>::New();
}
}
if (newNumPoints>0)
vtkNewCells->InsertNextCell(container);
}
}
if (vtkNewCells->GetNumberOfCells()<=0)
return nullptr;
vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::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<PlanarFigure*>(roi->GetData()) || dynamic_cast<PlanarFigureComposite*>(roi->GetData())) )
return nullptr;
std::vector<long> tmp = ExtractFiberIdSubset(roi, storage);
if (tmp.size()<=0)
return mitk::FiberBundle::New();
vtkSmartPointer<vtkPolyData> pTmp = GeneratePolyDataByIds(tmp);
return mitk::FiberBundle::New(pTmp);
}
std::vector<long> mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage)
{
std::vector<long> result;
if (roi==nullptr || roi->GetData()==nullptr)
return result;
mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast<mitk::PlanarFigureComposite*>(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<long>::iterator it;
for (unsigned int i=1; i<children->Size(); ++i)
{
std::vector<long> inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage);
std::vector<long> 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<long>::iterator it;
for (unsigned int i=1; i<children->Size(); ++i)
{
it = result.end();
std::vector<long> 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; i<this->GetNumFibers(); i++)
result.push_back(i);
std::vector<long>::iterator it;
for (unsigned int i=0; i<children->Size(); ++i)
{
std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
std::vector<long> 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<mitk::PlanarFigure*>(roi->GetData()) ) // actual extraction
{
if ( dynamic_cast<mitk::PlanarPolygon*>(roi->GetData()) )
{
mitk::PlanarFigure::Pointer planarPoly = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
//create vtkPolygon using controlpoints from planarFigure polygon
vtkSmartPointer<vtkPolygon> polygonVtk = vtkSmartPointer<vtkPolygon>::New();
for (unsigned int i=0; i<planarPoly->GetNumberOfControlPoints(); ++i)
{
itk::Point<double,3> 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; i<m_NumFibers; i++)
{
++disp ;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
for (int j=0; j<numPoints-1; j++)
{
// Inputs
double p1[3] = {0,0,0};
points->GetPoint(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<mitk::PlanarCircle*>(roi->GetData()) )
{
mitk::PlanarFigure::Pointer planarFigure = dynamic_cast<mitk::PlanarFigure*>(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; i<m_NumFibers; i++)
{
++disp ;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
for (int j=0; j<numPoints-1; j++)
{
// Inputs
double p1[3] = {0,0,0};
points->GetPoint(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<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::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<vtkFloatArray>::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; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int p = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
float length = 0;
for (int j=0; j<p-1; j++)
{
double p1[3];
points->GetPoint(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 (length<m_MinFiberLength)
m_MinFiberLength = length;
if (length>m_MaxFiberLength)
m_MaxFiberLength = length;
}
}
m_MeanFiberLength /= m_NumFibers;
std::vector< float > sortedLengths = m_FiberLengths;
std::sort(sortedLengths.begin(), sortedLengths.end());
for (int i=0; i<m_NumFibers; i++)
m_LengthStDev += (m_MeanFiberLength-sortedLengths.at(i))*(m_MeanFiberLength-sortedLengths.at(i));
if (m_NumFibers>1)
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; i<m_FiberWeights->GetSize(); i++)
+ for (int i=0; i<m_FiberWeights->GetNumberOfValues(); i++)
m_FiberWeights->SetValue(i, newWeight);
}
void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer<vtkFloatArray> 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; i<weights->GetSize(); i++)
+ for (int i=0; i<weights->GetNumberOfValues(); 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<vtkUnsignedCharArray> fiberColors)
{
for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++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<float, 3> 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<float, 3> 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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
for (int i=0; i<m_NumFibers; i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(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<vtkPolyData>::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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
for (int i=0; i<m_NumFibers; i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(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<vtkPolyData>::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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
for (int i=0; i<m_NumFibers; i++)
{
++disp ;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(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<vtkPolyData>::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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
for (int i=0; i<m_NumFibers; i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(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<vtkPolyData>::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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
for (int i=0; i<m_NumFibers; i++)
{
++disp;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
p[axis] = -p[axis];
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
vtkNewCells->InsertNextCell(container);
}
m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
m_FiberPolyData->SetPoints(vtkNewPoints);
m_FiberPolyData->SetLines(vtkNewCells);
this->SetFiberPolyData(m_FiberPolyData, true);
}
void mitk::FiberBundle::RemoveDir(vnl_vector_fixed<double,3> dir, double threshold)
{
dir.normalize();
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
++disp ;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
// calculate curvatures
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
bool discard = false;
for (int j=0; j<numPoints-1; j++)
{
double p1[3];
points->GetPoint(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; j<numPoints; j++)
{
double p1[3];
points->GetPoint(j, p1);
vtkIdType id = vtkNewPoints->InsertNextPoint(p1);
container->GetPointIds()->InsertNextId(id);
}
vtkNewCells->InsertNextCell(container);
}
}
m_FiberPolyData = vtkSmartPointer<vtkPolyData>::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<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
MITK_INFO << "Applying curvature threshold";
boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
++disp ;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
// calculate curvatures
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints-2; j++)
{
double p1[3];
points->GetPoint(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 && r<minRadius)
break;
if (r<minRadius)
{
j += 2;
vtkNewCells->InsertNextCell(container);
container = vtkSmartPointer<vtkPolyLine>::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<vtkPolyData>::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 || lengthInMM<m_MinFiberLength)
{
MITK_INFO << "No fibers shorter than " << lengthInMM << " mm found!";
return true;
}
if (lengthInMM>m_MaxFiberLength) // can't remove all fibers
{
MITK_WARN << "Process aborted. No fibers would be left!";
return false;
}
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
float min = m_MaxFiberLength;
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
if (m_FiberLengths.at(i)>=lengthInMM)
{
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
vtkNewCells->InsertNextCell(container);
if (m_FiberLengths.at(i)<min)
min = m_FiberLengths.at(i);
}
}
if (vtkNewCells->GetNumberOfCells()<=0)
return false;
m_FiberPolyData = vtkSmartPointer<vtkPolyData>::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<m_MinFiberLength) // can't remove all fibers
return false;
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
MITK_INFO << "Removing long fibers";
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
if (m_FiberLengths.at(i)<=lengthInMM)
{
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
vtkNewCells->InsertNextCell(container);
}
}
if (vtkNewCells->GetNumberOfCells()<=0)
return false;
m_FiberPolyData = vtkSmartPointer<vtkPolyData>::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<vtkPoints> vtkSmoothPoints = vtkSmartPointer<vtkPoints>::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<vtkCellArray> vtkSmoothCells = vtkSmartPointer<vtkCellArray>::New(); //cellcontainer for smoothed lines
vtkIdType pointHelperCnt = 0;
MITK_INFO << "Smoothing fibers";
vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
newFiberWeights->SetName("FIBER_WEIGHTS");
newFiberWeights->SetNumberOfValues(m_NumFibers);
boost::progress_display disp(m_NumFibers);
#pragma omp parallel for
for (int i=0; i<m_NumFibers; i++)
{
vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::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; j<numPoints; j++)
newPoints->InsertNextPoint(points->GetPoint(j));
}
int sampling = std::ceil(length/pointDistance);
vtkSmartPointer<vtkKochanekSpline> xSpline = vtkSmartPointer<vtkKochanekSpline>::New();
vtkSmartPointer<vtkKochanekSpline> ySpline = vtkSmartPointer<vtkKochanekSpline>::New();
vtkSmartPointer<vtkKochanekSpline> zSpline = vtkSmartPointer<vtkKochanekSpline>::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<vtkParametricSpline> spline = vtkSmartPointer<vtkParametricSpline>::New();
spline->SetXSpline(xSpline);
spline->SetYSpline(ySpline);
spline->SetZSpline(zSpline);
spline->SetPoints(newPoints);
vtkSmartPointer<vtkParametricFunctionSource> functionSource = vtkSmartPointer<vtkParametricFunctionSource>::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<vtkPolyLine> smoothLine = vtkSmartPointer<vtkPolyLine>::New();
smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints());
#pragma omp critical
{
for (int j=0; j<smoothLine->GetNumberOfPoints(); 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<vtkPolyData>::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; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
points += cell->GetNumberOfPoints();
}
return points;
}
void mitk::FiberBundle::Compress(float error)
{
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
MITK_INFO << "Compressing fibers";
unsigned long numRemovedPoints = 0;
boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
newFiberWeights->SetName("FIBER_WEIGHTS");
newFiberWeights->SetNumberOfValues(m_NumFibers);
#pragma omp parallel for
for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); 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; j<numPoints; j++)
{
double cand[3];
points->GetPoint(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<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
int remCounter = 0;
bool pointFound = true;
while (pointFound)
{
pointFound = false;
double minError = error;
int removeIndex = -1;
for (unsigned int j=0; j<vertices.size(); j++)
{
if (removedPoints[j]==0)
{
vnl_vector_fixed< double, 3 > 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<numPoints; k++)
if (removedPoints[k]<=0)
{
succ = vertices.at(k);
validS = k;
break;
}
if (validP>=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 (hc<minError)
{
removeIndex = j;
minError = hc;
pointFound = true;
}
}
}
}
if (pointFound)
{
removedPoints[removeIndex] = 1;
remCounter++;
}
}
for (int j=0; j<numPoints; j++)
{
if (removedPoints[j]<=0)
{
#pragma omp critical
{
vtkIdType id = vtkNewPoints->InsertNextPoint(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<vtkPolyData>::New();
m_FiberPolyData->SetPoints(vtkNewPoints);
m_FiberPolyData->SetLines(vtkNewCells);
this->SetFiberPolyData(m_FiberPolyData, true);
}
}
void mitk::FiberBundle::ResampleLinear(double pointDistance)
{
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
MITK_INFO << "Resampling fibers (linear)";
boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
vtkSmartPointer<vtkFloatArray> newFiberWeights = vtkSmartPointer<vtkFloatArray>::New();
newFiberWeights->SetName("FIBER_WEIGHTS");
newFiberWeights->SetNumberOfValues(m_NumFibers);
#pragma omp parallel for
for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); 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; j<numPoints; j++)
{
double cand[3];
points->GetPoint(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<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::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<vertices.size(); j++)
{
vnl_vector_fixed< double, 3 > 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<vtkPolyData>::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; i<m_NumFibers; i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(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; j<numPoints; j++)
{
double* p1 = points->GetPoint(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 <mitkBaseData.h>
#include <MitkFiberTrackingExports.h>
#include <mitkImage.h>
#include <mitkDataStorage.h>
#include <mitkPlanarFigure.h>
#include <mitkPixelTypeTraits.h>
#include <mitkPlanarFigureComposite.h>
//includes storing fiberdata
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkPoints.h>
#include <vtkDataSet.h>
#include <vtkTransform.h>
#include <vtkFloatArray.h>
namespace mitk {
/**
* \brief Base Class for Fiber Bundles; */
class MITKFIBERTRACKING_EXPORT FiberBundle : public BaseData
{
public:
typedef itk::Image<unsigned char, 3> 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<vtkPolyData>) // 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 <typename TPixel>
void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity, bool normalize);
void ColorFibersByOrientation();
void SetFiberOpacity(vtkDoubleArray *FAValArray);
void ResetFiberOpacity();
void SetFiberColors(vtkSmartPointer<vtkUnsignedCharArray> fiberColors);
void SetFiberColors(float r, float g, float b, float alpha=255);
vtkSmartPointer<vtkUnsignedCharArray> 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<double,3> dir, double threshold);
itk::Point<float, 3> 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<long> 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<vtkPolyData> GeneratePolyDataByIds( std::vector<long> ); // TODO: make protected
void GenerateFiberIds(); // TODO: make protected
// get/set data
vtkSmartPointer<vtkFloatArray> 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<vtkFloatArray> weights);
void SetFiberPolyData(vtkSmartPointer<vtkPolyData>, bool updateGeometry = true);
vtkSmartPointer<vtkPolyData> 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<float, 3> 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<vtkPolyData> m_FiberPolyData;
// contains fiber ids
vtkSmartPointer<vtkDataSet> m_FiberIdDataSet;
int m_NumFibers;
vtkSmartPointer<vtkUnsignedCharArray> m_FiberColors;
vtkSmartPointer<vtkFloatArray> 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 <mitkImageCast.h>
#include <mitkITKImageImport.h>
#include <mitkProperties.h>
#include <mitkImage.h>
#include <mitkIOUtil.h>
#include <mitkFiberBundle.h>
#include <mitkFiberfoxParameters.h>
#include "mitkCommandLineParser.h"
#include <itkTractsToDWIImageFilter.h>
#include <boost/lexical_cast.hpp>
#include <mitkPreferenceListReaderOptionsFunctor.h>
#include <itksys/SystemTools.hxx>
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<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
{
return EXIT_FAILURE;
}
string outName = us::any_cast<string>(parsedArgs["out"]);
string paramName = us::any_cast<string>(parsedArgs["parameters"]);
string input="";
if (parsedArgs.count("input"))
{
input = us::any_cast<string>(parsedArgs["input"]);
}
bool verbose = false;
if (parsedArgs.count("verbose"))
{
verbose = us::any_cast<bool>(parsedArgs["verbose"]);
}
FiberfoxParameters<double> 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<int>::max() )
{
outName = NameTest + "_" + boost::lexical_cast<std::string>(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<mitk::FiberBundle*>(inputData.GetPointer()) ) // simulate dataset from fibers
{
tractsToDwiFilter->SetFiberBundle(dynamic_cast<mitk::FiberBundle*>(inputData.GetPointer()));
}
else if ( dynamic_cast<mitk::Image*>(inputData.GetPointer()) ) // add artifacts to existing image
{
typedef itk::VectorImage< short, 3 > ItkDwiType;
mitk::Image::Pointer diffImg = dynamic_cast<mitk::Image*>(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<mitk::FloatProperty*>
(diffImg->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )
->GetValue();
parameters.m_SignalGen.SetGradienDirections( static_cast<mitk::GradientDirectionsProperty*>
( 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; k<volumeFractions.size(); k++)
{
mitk::Image::Pointer image = mitk::Image::New();
image->InitializeByItk(volumeFractions.at(k).GetPointer());
image->SetVolume(volumeFractions.at(k)->GetBufferPointer());
mitk::IOUtil::Save(image, outName+"_Compartment"+boost::lexical_cast<string>(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 <mitkBaseData.h>
+#include <mitkImageCast.h>
+#include <mitkImageToItk.h>
+#include <metaCommand.h>
+#include <mitkCommandLineParser.h>
+#include <usAny.h>
+#include <mitkIOUtil.h>
+#include <boost/lexical_cast.hpp>
+#include <itksys/SystemTools.hxx>
+#include <itkDirectory.h>
+#include <mitkFiberBundle.h>
+#include <mitkPreferenceListReaderOptionsFunctor.h>
+#include <mitkDiffusionPropertyHelper.h>
+#include <vnl/vnl_linear_system.h>
+#include <Eigen/Dense>
+#include <mitkStickModel.h>
+#include <mitkBallModel.h>
+#include <vigra/regression.hxx>
+#include <itkImageFileWriter.h>
+#include <itkImageDuplicator.h>
+#include <itkMersenneTwisterRandomVariateGenerator.h>
+
+using namespace std;
+typedef itksys::SystemTools ist;
+typedef itk::Image<unsigned char, 3> ItkUcharImgType;
+typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType;
+typedef mitk::DiffusionPropertyHelper DPH;
+typedef itk::Point<float, 3> PointType;
+typedef mitk::StickModel<> ModelType;
+typedef mitk::BallModel<> BallModelType;
+
+std::vector<float> 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<vtkPolyData> polydata = inputTractogram->GetFiberPolyData();
+ for (int i=0; i<inputTractogram->GetNumFibers(); ++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; j<numPoints-1; ++j)
+ {
+ double* p1 = points->GetPoint(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; k<num_gradients; ++k)
+ {
+ b(num_gradients*linear_index + k, 0) = (double)measured_signal[k];
+ A(num_gradients*linear_index + k, i) = A(num_gradients*linear_index + k, i) + (double)model_signal[k];
+
+ temp_img1->SetPixel(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<float> weights;
+ for (unsigned int i=0; i<num_unknowns; ++i)
+ weights.push_back(x[i]);
+ return weights;
+}
+
+std::vector<float> SolveEvo(mitk::FiberBundle* inputTractogram, mitk::Image::Pointer inputImage, ModelType signalModel, BallModelType ballModel, int num_iterations=1000, double step=0.1)
+{
+ std::vector<float> 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<vtkPolyData> 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; i<inputTractogram->GetNumFibers(); ++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; j<numPoints-1; ++j)
+ {
+ double* p1 = points->GetPoint(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<num_iterations; ++i)
+ {
+ MITK_INFO << "Iteration " << i;
+ std::random_shuffle(fiber_indices.begin(), fiber_indices.end());
+
+ double T = T_start * exp(alpha*(double)i/num_iterations);
+
+ int accepted = 0;
+ for (auto f : fiber_indices)
+ {
+ std::vector< ModelType::PixelType > 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]<step)
+ add = 1;
+
+ int c = 0;
+ for (auto idx : image_indices)
+ {
+ DPH::ImageType::PixelType mVal = itkImage->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<num_gradients; ++g)
+ {
+// ModelType::GradientType d = signalModel.GetGradientDirection(g);
+// if (d.GetNorm()<0.0001)
+// continue;
+
+ double delta = (double)mVal[g] - sVal[g];
+
+ E_old += fabs(delta)/num_gradients;
+
+ if (add==1)
+ {
+ delta = (double)mVal[g] - (sVal[g] + dVal[g]);
+ E_new += fabs(delta)/num_gradients;
+ }
+ else
+ {
+ delta = (double)mVal[g] - (sVal[g] - dVal[g]);
+ E_new += fabs(delta)/num_gradients;
+ }
+ }
+ ++c;
+ }
+
+ E_old /= image_indices.size();
+ E_new /= image_indices.size();
+ double R = exp( (E_old-E_new)/T );
+
+// if (E_new < E_old)
+ float p = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
+ if (p<R)
+ {
+ accepted++;
+ c = 0;
+ for (auto idx : image_indices)
+ {
+ ModelType::PixelType sVal = simulatedImage->GetPixel(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<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+ if (parsedArgs.size()==0)
+ return EXIT_FAILURE;
+
+ string fibFile = us::any_cast<string>(parsedArgs["input_tractogram"]);
+ string dwiFile = us::any_cast<string>(parsedArgs["input_dMRI"]);
+ string outRoot = us::any_cast<string>(parsedArgs["out"]);
+
+ int it = 1000;
+ if (parsedArgs.count("it"))
+ it = us::any_cast<int>(parsedArgs["it"]);
+
+ float step = 0.1;
+ if (parsedArgs.count("step"))
+ step = us::any_cast<float>(parsedArgs["step"]);
+
+
+
+ typedef DPH::GradientDirectionsContainerType GradientContainerType;
+
+ try
+ {
+ mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fibFile)[0].GetPointer());
+
+ mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "Fiberbundles"}, {});
+ mitk::Image::Pointer inputImage = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(dwiFile, &functor)[0].GetPointer());
+ if (!DPH::IsDiffusionWeightedImage(inputImage))
+ return EXIT_FAILURE;
+
+ // resample fibers
+ float minSpacing = 1;
+ if(inputImage->GetGeometry()->GetSpacing()[0]<inputImage->GetGeometry()->GetSpacing()[1] && inputImage->GetGeometry()->GetSpacing()[0]<inputImage->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<mitk::GradientDirectionsProperty*>( 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<float> weights = SolveLinear(inputTractogram, inputImage, signalModel, ballModel);
+ std::vector<float> weights = SolveEvo(inputTractogram, inputImage, signalModel, ballModel, it, step);
+
+ for (unsigned int i=0; i<weights.size(); ++i)
+ {
+ MITK_INFO << weights.at(i);
+ inputTractogram->SetFiberWeight(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 <mitkBaseData.h>
+#include <mitkImageCast.h>
+#include <mitkImageToItk.h>
+#include <metaCommand.h>
+#include <mitkCommandLineParser.h>
+#include <usAny.h>
+#include <mitkIOUtil.h>
+#include <boost/lexical_cast.hpp>
+#include <itksys/SystemTools.hxx>
+#include <itkDirectory.h>
+#include <mitkFiberBundle.h>
+#include <vtkTransformPolyDataFilter.h>
+
+using namespace std;
+typedef itksys::SystemTools ist;
+typedef itk::Image<unsigned char, 3> 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::Image*>(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<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::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<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+ if (parsedArgs.size()==0)
+ return EXIT_FAILURE;
+
+ string fibFile = us::any_cast<string>(parsedArgs["input"]);
+ string gray_matter_mask = us::any_cast<string>(parsedArgs["gray_matter_mask"]);
+ string reference_mask_folder = us::any_cast<string>(parsedArgs["reference_mask_folder"]);
+ string outRoot = us::any_cast<string>(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::FiberBundle*>(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]<std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] && std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]<std::get<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 <berryISelectionService.h>
#include <berryIWorkbenchPart.h>
#include <berryIWorkbenchWindow.h>
// Qmitk
#include "QmitkFiberProcessingView.h"
// Qt
#include <QMessageBox>
// MITK
#include <mitkNodePredicateProperty.h>
#include <mitkImageCast.h>
#include <mitkPointSet.h>
#include <mitkPlanarCircle.h>
#include <mitkPlanarPolygon.h>
#include <mitkPlanarRectangle.h>
#include <mitkPlanarFigureInteractor.h>
#include <mitkImageAccessByItk.h>
#include <mitkDataNodeObject.h>
#include <mitkTensorImage.h>
#include "usModuleRegistry.h"
#include <itkFiberCurvatureFilter.h>
#include "mitkNodePredicateDataType.h"
#include <mitkNodePredicateProperty.h>
#include <mitkNodePredicateAnd.h>
#include <mitkNodePredicateNot.h>
#include <mitkNodePredicateOr.h>
// ITK
#include <itkResampleImageFilter.h>
#include <itkGaussianInterpolateImageFunction.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkTractsToFiberEndingsImageFilter.h>
#include <itkTractDensityImageFilter.h>
#include <itkImageRegion.h>
#include <itkTractsToRgbaImageFilter.h>
#define _USE_MATH_DEFINES
#include <math.h>
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<mitk::Image>::Pointer isMitkImage = mitk::TNodePredicateDataType<mitk::Image>::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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(node->GetData());
vnl_vector_fixed<double,3> 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<mitk::Image*>(m_MaskImageNode->GetData());
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(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<mitk::Image*>(m_MaskImageNode->GetData());
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(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<mitk::FiberBundle*>(m_SelectedFB.front()->GetData());
geometry = fib->GetGeometry();
}
else if (m_SelectedImage)
geometry = m_SelectedImage->GetGeometry();
else
return;
itk::Vector<double,3> 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<double, 3, 3> 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; i<m_SelectedPF.size(); i++)
{
name += "+";
name += m_SelectedPF.at(i)->GetName();
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<mitk::PlanarFigure*>(node->GetData()))
{
m_PlanarFigure = dynamic_cast<mitk::PlanarFigure*>(node->GetData());
AccessFixedDimensionByItk_2(
image,
InternalReorientImagePlane, 3,
m_PlanarFigure->GetGeometry(), -1);
AccessFixedDimensionByItk_2(
m_InternalImage,
InternalCalculateMaskFromPlanarFigure,
3, 2, node->GetName() );
}
else if (dynamic_cast<mitk::PlanarFigureComposite*>(node->GetData()))
{
DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node);
for (unsigned int i=0; i<children->Size(); 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<ImageType, FloatImageType, double> ResamplerType;
typename ResamplerType::Pointer resampler = ResamplerType::New();
mitk::PlaneGeometry* planegeo = dynamic_cast<mitk::PlaneGeometry*>(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; c<matrix.ColumnDimensions; c++)
{
double sum = 0;
for(unsigned int r=0; r<matrix.RowDimensions; r++)
sum += matrix(r,c)*matrix(r,c);
for(unsigned int r=0; r<matrix.RowDimensions; r++)
direction(r,c) = matrix(r,c)/sqrt(sum);
}
resampler->SetOutputDirection( 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<ImageType, double> GaussianInterpolatorType;
typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New();
interpolator->SetInputImage( image );
interpolator->SetParameters( sigma, alpha );
resampler->SetInterpolator( interpolator );
}
else
{
typedef typename itk::LinearInterpolateImageFunction<ImageType, double> 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<itkUCharImageType>
itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion());
itk::ImageRegionIterator<ImageType>
itimage(image, image->GetLargestPossibleRegion());
itmask.GoToBegin();
itimage.GoToBegin();
typename ImageType::SizeType lowersize = {{itk::NumericTraits<unsigned long>::max(),itk::NumericTraits<unsigned long>::max(),itk::NumericTraits<unsigned long>::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<itkUCharImageType> 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("<font color='red'>mandatory</font>");
m_Controls->m_PfLabel->setText("<font color='grey'>needed for extraction</font>");
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<mitk::DataNode::Pointer> nodes;
OnSelectionChanged(nullPart, nodes);
}
void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* )
{
if (!m_Controls->m_InteractiveBox->isChecked())
{
berry::IWorkbenchPart::Pointer nullPart;
QList<mitk::DataNode::Pointer> 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<mitk::PlanarFigure*>(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<mitk::PlanarFigure*>(node->GetData());
if (figure!=nullptr)
figure->RemoveAllObservers();
}
}
void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList<mitk::DataNode::Pointer>& nodes)
{
RemoveObservers();
//reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection
std::vector<mitk::DataNode::Pointer> 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<mitk::FiberBundle*>(node->GetData()) )
m_SelectedFB.push_back(node);
else if (dynamic_cast<mitk::PlanarPolygon*>(node->GetData()) || dynamic_cast<mitk::PlanarFigureComposite*>(node->GetData()) || dynamic_cast<mitk::PlanarCircle*>(node->GetData()))
m_SelectedPF.push_back(node);
else if (dynamic_cast<mitk::Image*>(node->GetData()))
{
m_SelectedImage = dynamic_cast<mitk::Image*>(node->GetData());
bool isBinary = false;
node->GetPropertyValue<bool>("binary", isBinary);
if (isBinary)
m_MaskImageNode = node;
}
else if (dynamic_cast<mitk::Surface*>(node->GetData()))
m_SelectedSurfaces.push_back(dynamic_cast<mitk::Surface*>(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<unsigned int, mitk::DataNode::Pointer>::ConstPointer nodes = this->GetDataStorage()->GetAll();
for (unsigned int i=0; i<nodes->Size(); i++)
if (dynamic_cast<mitk::FiberBundle*>(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<unsigned int, mitk::DataNode::Pointer>::ConstPointer nodes = this->GetDataStorage()->GetAll();
for (unsigned int i=0; i<nodes->Size(); i++)
if (dynamic_cast<mitk::PlanarPolygon*>(nodes->at(i)->GetData()) || dynamic_cast<mitk::PlanarFigureComposite*>(nodes->at(i)->GetData()) || dynamic_cast<mitk::PlanarCircle*>(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<mitk::PlanarFigureInteractor*>(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<mitk::DataNode::Pointer> fiberBundles = m_SelectedFB;
mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0);
for (unsigned int i=0; i<fiberBundles.size(); i++)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(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<mitk::DataNode::Pointer> 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<PlanarFigure*>(child->GetData()))
{
mitk::DataNode::Pointer newChild;
newChild = mitk::DataNode::New();
newChild->SetData(dynamic_cast<PlanarFigure*>(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<PlanarFigureComposite*>(child->GetData()))
{
mitk::DataNode::Pointer newChild;
newChild = mitk::DataNode::New();
newChild->SetData(dynamic_cast<PlanarFigureComposite*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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; i<m_SelectedFB.size(); i++)
{
newBundle = newBundle->AddBundle(dynamic_cast<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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; i<m_SelectedFB.size(); i++)
{
newBundle = newBundle->SubtractBundle(dynamic_cast<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(node->GetData());
fib->ColorFibersByScalarMap(dynamic_cast<mitk::Image*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(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<mitk::FiberBundle*>(node->GetData());
if (m_SelectedImage.IsNotNull())
fib->SetReferenceGeometry(m_SelectedImage->GetGeometry());
fib->MirrorFibers(axis);
}
for (auto surf : m_SelectedSurfaces)
{
vtkSmartPointer<vtkPolyData> poly = surf->GetVtkPolyData();
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
for (int i=0; i<poly->GetNumberOfPoints(); i++)
{
double* point = poly->GetPoint(i);
point[axis] *= -1;
vtkNewPoints->InsertNextPoint(point);
}
poly->SetPoints(vtkNewPoints);
surf->CalculateBoundingBox();
}
if (auto renderWindowPart = this->GetRenderWindowPart())
{
renderWindowPart->RequestUpdate();
}
}