diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
index 9b0a9ed25e..d9a719ec3d 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp
@@ -1,427 +1,437 @@
#include "itkFitFibersToImageFilter.h"
#include <boost/progress.hpp>
namespace itk{
FitFibersToImageFilter::FitFibersToImageFilter()
- : m_FitIndividualFibers(true)
+ : m_PeakImage(nullptr)
+ , m_MaskImage(nullptr)
+ , m_FitIndividualFibers(true)
, m_GradientTolerance(1e-5)
, m_Lambda(0.1)
, m_MaxIterations(20)
, m_FiberSampling(10)
, m_Coverage(0)
, m_Overshoot(0)
, m_FilterOutliers(true)
, m_MeanWeight(1.0)
, m_MedianWeight(1.0)
, m_MinWeight(1.0)
, m_MaxWeight(1.0)
, m_Verbose(true)
, m_DeepCopy(true)
+ , m_NumUnknowns(0)
+ , m_NumResiduals(0)
+ , m_NumCoveredDirections(0)
{
this->SetNumberOfRequiredOutputs(3);
}
FitFibersToImageFilter::~FitFibersToImageFilter()
{
}
void FitFibersToImageFilter::GenerateData()
{
int sz_x = m_PeakImage->GetLargestPossibleRegion().GetSize(0);
int sz_y = m_PeakImage->GetLargestPossibleRegion().GetSize(1);
int sz_z = m_PeakImage->GetLargestPossibleRegion().GetSize(2);
int sz_peaks = m_PeakImage->GetLargestPossibleRegion().GetSize(3)/3 + 1; // +1 for zero - peak
int num_voxels = sz_x*sz_y*sz_z;
float minSpacing = 1;
if(m_PeakImage->GetSpacing()[0]<m_PeakImage->GetSpacing()[1] && m_PeakImage->GetSpacing()[0]<m_PeakImage->GetSpacing()[2])
minSpacing = m_PeakImage->GetSpacing()[0];
else if (m_PeakImage->GetSpacing()[1] < m_PeakImage->GetSpacing()[2])
minSpacing = m_PeakImage->GetSpacing()[1];
else
minSpacing = m_PeakImage->GetSpacing()[2];
- unsigned int num_unknowns = m_Tractograms.size();
+ m_NumUnknowns = m_Tractograms.size();
if (m_FitIndividualFibers)
{
- num_unknowns = 0;
+ m_NumUnknowns = 0;
for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
- num_unknowns += m_Tractograms.at(bundle)->GetNumFibers();
+ m_NumUnknowns += m_Tractograms.at(bundle)->GetNumFibers();
}
for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
{
if (m_DeepCopy)
m_Tractograms[bundle] = m_Tractograms.at(bundle)->GetDeepCopy();
m_Tractograms.at(bundle)->ResampleLinear(minSpacing/m_FiberSampling);
}
- unsigned int number_of_residuals = num_voxels * sz_peaks;
+ m_NumResiduals = num_voxels * sz_peaks;
- MITK_INFO << "Num. unknowns: " << num_unknowns;
- MITK_INFO << "Num. residuals: " << number_of_residuals;
+ MITK_INFO << "Num. unknowns: " << m_NumUnknowns;
+ MITK_INFO << "Num. residuals: " << m_NumResiduals;
MITK_INFO << "Creating system ...";
vnl_sparse_matrix<double> A;
vnl_vector<double> b;
- A.set_size(number_of_residuals, num_unknowns);
- b.set_size(number_of_residuals); b.fill(0.0);
+ A.set_size(m_NumResiduals, m_NumUnknowns);
+ b.set_size(m_NumResiduals); b.fill(0.0);
double TD = 0;
double FD = 0;
- unsigned int dir_count = 0;
+ m_NumCoveredDirections = 0;
unsigned int fiber_count = 0;
for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
{
vtkSmartPointer<vtkPolyData> polydata = m_Tractograms.at(bundle)->GetFiberPolyData();
for (int i=0; i<m_Tractograms.at(bundle)->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);
PointType4 p;
p[0]=p1[0];
p[1]=p1[1];
p[2]=p1[2];
p[3]=0;
itk::Index<4> idx4;
m_PeakImage->TransformPhysicalPointToIndex(p, idx4);
- if (!m_PeakImage->GetLargestPossibleRegion().IsInside(idx4))
+ itk::Index<3> idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2];
+ if (!m_PeakImage->GetLargestPossibleRegion().IsInside(idx4) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0))
continue;
double* p2 = points->GetPoint(j+1);
vnl_vector_fixed<float,3> fiber_dir;
fiber_dir[0] = p[0]-p2[0];
fiber_dir[1] = p[1]-p2[1];
fiber_dir[2] = p[2]-p2[2];
fiber_dir.normalize();
double w = 1;
int peak_id = sz_peaks-1;
vnl_vector_fixed<float,3> odf_peak = GetClosestPeak(idx4, m_PeakImage, fiber_dir, peak_id, w);
float peak_mag = odf_peak.magnitude();
int x = idx4[0];
int y = idx4[1];
int z = idx4[2];
unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z + sz_x*sz_y*sz_z*peak_id;
if (b[linear_index] == 0 && peak_id<3)
{
- dir_count++;
+ m_NumCoveredDirections++;
FD += peak_mag;
}
TD += w;
if (m_FitIndividualFibers)
{
b[linear_index] = peak_mag;
A.put(linear_index, fiber_count, A.get(linear_index, fiber_count) + w);
}
else
{
b[linear_index] = peak_mag;
A.put(linear_index, bundle, A.get(linear_index, bundle) + w);
}
}
++fiber_count;
}
}
- TD /= (dir_count*fiber_count);
- FD /= dir_count;
+ TD /= (m_NumCoveredDirections*fiber_count);
+ FD /= m_NumCoveredDirections;
A /= TD;
b *= 100.0/FD; // times 100 because we want to avoid too small values for computational reasons
double init_lambda = fiber_count; // initialization for lambda estimation
itk::TimeProbe clock;
clock.Start();
- VnlCostFunction cost(num_unknowns);
+ VnlCostFunction cost(m_NumUnknowns);
cost.SetProblem(A, b, init_lambda);
- m_Weights.set_size(num_unknowns); // m_Weights.fill( TD/100.0 * FD/2.0 );
+ m_Weights.set_size(m_NumUnknowns); // m_Weights.fill( TD/100.0 * FD/2.0 );
m_Weights.fill( 0.0 );
vnl_lbfgsb minimizer(cost);
- vnl_vector<double> l; l.set_size(num_unknowns); l.fill(0);
- vnl_vector<long> bound_selection; bound_selection.set_size(num_unknowns); bound_selection.fill(1);
+ vnl_vector<double> l; l.set_size(m_NumUnknowns); l.fill(0);
+ vnl_vector<long> bound_selection; bound_selection.set_size(m_NumUnknowns); bound_selection.fill(1);
minimizer.set_bound_selection(bound_selection);
minimizer.set_lower_bound(l);
minimizer.set_projected_gradient_tolerance(m_GradientTolerance);
MITK_INFO << "Estimating regularization";
minimizer.set_trace(false);
minimizer.set_max_function_evals(2);
minimizer.minimize(m_Weights);
- vnl_vector<double> dx; dx.set_size(num_unknowns); dx.fill(0.0);
+ vnl_vector<double> dx; dx.set_size(m_NumUnknowns); dx.fill(0.0);
cost.grad_regu_localMSE(m_Weights, dx);
double r = dx.magnitude()/m_Weights.magnitude();
cost.m_Lambda *= m_Lambda*55.0/r;
MITK_INFO << r << " - " << m_Lambda*55.0/r;
if (cost.m_Lambda>10e7)
{
MITK_INFO << "Regularization estimation failed. Using default value.";
cost.m_Lambda = fiber_count;
}
MITK_INFO << "Using regularization factor of " << cost.m_Lambda << " (λ: " << m_Lambda << ")";
MITK_INFO << "Fitting fibers";
minimizer.set_trace(m_Verbose);
minimizer.set_max_function_evals(m_MaxIterations);
minimizer.minimize(m_Weights);
std::vector< double > weights;
if (m_FilterOutliers)
{
for (auto w : m_Weights)
weights.push_back(w);
std::sort(weights.begin(), weights.end());
- MITK_INFO << "Setting upper weight bound to " << weights.at(num_unknowns*0.99);
- vnl_vector<double> u; u.set_size(num_unknowns); u.fill(weights.at(num_unknowns*0.99));
+ MITK_INFO << "Setting upper weight bound to " << weights.at(m_NumUnknowns*0.99);
+ vnl_vector<double> u; u.set_size(m_NumUnknowns); u.fill(weights.at(m_NumUnknowns*0.99));
minimizer.set_upper_bound(u);
bound_selection.fill(2);
minimizer.set_bound_selection(bound_selection);
minimizer.minimize(m_Weights);
weights.clear();
}
for (auto w : m_Weights)
weights.push_back(w);
std::sort(weights.begin(), weights.end());
m_MeanWeight = m_Weights.mean();
- m_MedianWeight = weights.at(num_unknowns*0.5);
+ m_MedianWeight = weights.at(m_NumUnknowns*0.5);
m_MinWeight = weights.at(0);
- m_MaxWeight = weights.at(num_unknowns-1);
+ m_MaxWeight = weights.at(m_NumUnknowns-1);
MITK_INFO << "*************************";
MITK_INFO << "Weight statistics";
MITK_INFO << "Mean: " << m_MeanWeight;
MITK_INFO << "Median: " << m_MedianWeight;
- MITK_INFO << "75% quantile: " << weights.at(num_unknowns*0.75);
- MITK_INFO << "95% quantile: " << weights.at(num_unknowns*0.95);
- MITK_INFO << "99% quantile: " << weights.at(num_unknowns*0.99);
+ MITK_INFO << "75% quantile: " << weights.at(m_NumUnknowns*0.75);
+ MITK_INFO << "95% quantile: " << weights.at(m_NumUnknowns*0.95);
+ MITK_INFO << "99% quantile: " << weights.at(m_NumUnknowns*0.99);
MITK_INFO << "Min: " << m_MinWeight;
MITK_INFO << "Max: " << m_MaxWeight;
MITK_INFO << "*************************";
MITK_INFO << "NumEvals: " << minimizer.get_num_evaluations();
MITK_INFO << "NumIterations: " << minimizer.get_num_iterations();
MITK_INFO << "Residual cost: " << minimizer.get_end_error();
MITK_INFO << "Final RMS: " << cost.S->get_rms_error(m_Weights);
clock.Stop();
int h = clock.GetTotal()/3600;
int m = ((int)clock.GetTotal()%3600)/60;
int s = (int)clock.GetTotal()%60;
MITK_INFO << "Optimization took " << h << "h, " << m << "m and " << s << "s";
// transform back for peak image creation
A *= FD/100.0;
b *= FD/100.0;
MITK_INFO << "Weighting fibers";
if (m_FitIndividualFibers)
{
unsigned int fiber_count = 0;
for (unsigned int bundle=0; bundle<m_Tractograms.size(); bundle++)
{
for (int i=0; i<m_Tractograms.at(bundle)->GetNumFibers(); i++)
{
m_Tractograms.at(bundle)->SetFiberWeight(i, m_Weights[fiber_count]);
++fiber_count;
}
m_Tractograms.at(bundle)->Compress(0.1);
m_Tractograms.at(bundle)->ColorFibersByFiberWeights(false, true);
}
}
else
{
for (unsigned int i=0; i<m_Tractograms.size(); ++i)
{
m_Tractograms.at(i)->SetFiberWeights(m_Weights[i]);
m_Tractograms.at(i)->Compress(0.1);
m_Tractograms.at(i)->ColorFibersByFiberWeights(false, true);
}
}
MITK_INFO << "Generating output images ...";
itk::ImageDuplicator< PeakImgType >::Pointer duplicator = itk::ImageDuplicator< PeakImgType >::New();
duplicator->SetInputImage(m_PeakImage);
duplicator->Update();
m_UnderexplainedImage = duplicator->GetOutput();
m_UnderexplainedImage->FillBuffer(0.0);
duplicator->SetInputImage(m_UnderexplainedImage);
duplicator->Update();
m_OverexplainedImage = duplicator->GetOutput();
m_OverexplainedImage->FillBuffer(0.0);
duplicator->SetInputImage(m_OverexplainedImage);
duplicator->Update();
m_ResidualImage = duplicator->GetOutput();
m_ResidualImage->FillBuffer(0.0);
duplicator->SetInputImage(m_ResidualImage);
duplicator->Update();
m_FittedImage = duplicator->GetOutput();
m_FittedImage->FillBuffer(0.0);
vnl_vector<double> fitted_b; fitted_b.set_size(b.size());
cost.S->multiply(m_Weights, fitted_b);
for (unsigned int r=0; r<b.size(); r++)
{
- itk::Index<4> idx;
+ itk::Index<4> idx4;
unsigned int linear_index = r;
- idx[0] = linear_index % sz_x; linear_index /= sz_x;
- idx[1] = linear_index % sz_y; linear_index /= sz_y;
- idx[2] = linear_index % sz_z; linear_index /= sz_z;
+ idx4[0] = linear_index % sz_x; linear_index /= sz_x;
+ idx4[1] = linear_index % sz_y; linear_index /= sz_y;
+ idx4[2] = linear_index % sz_z; linear_index /= sz_z;
int peak_id = linear_index % sz_peaks;
if (peak_id<sz_peaks-1)
{
vnl_vector_fixed<float,3> peak_dir;
- idx[3] = peak_id*3;
- peak_dir[0] = m_PeakImage->GetPixel(idx);
- idx[3] += 1;
- peak_dir[1] = m_PeakImage->GetPixel(idx);
- idx[3] += 1;
- peak_dir[2] = m_PeakImage->GetPixel(idx);
+ idx4[3] = peak_id*3;
+ peak_dir[0] = m_PeakImage->GetPixel(idx4);
+ idx4[3] += 1;
+ peak_dir[1] = m_PeakImage->GetPixel(idx4);
+ idx4[3] += 1;
+ peak_dir[2] = m_PeakImage->GetPixel(idx4);
peak_dir.normalize();
peak_dir *= fitted_b[r];
- idx[3] = peak_id*3;
- m_FittedImage->SetPixel(idx, peak_dir[0]);
+ idx4[3] = peak_id*3;
+ m_FittedImage->SetPixel(idx4, peak_dir[0]);
- idx[3] += 1;
- m_FittedImage->SetPixel(idx, peak_dir[1]);
+ idx4[3] += 1;
+ m_FittedImage->SetPixel(idx4, peak_dir[1]);
- idx[3] += 1;
- m_FittedImage->SetPixel(idx, peak_dir[2]);
+ idx4[3] += 1;
+ m_FittedImage->SetPixel(idx4, peak_dir[2]);
}
}
FD = 0;
m_Coverage = 0;
m_Overshoot = 0;
- itk::Index<4> idx;
- for (idx[0]=0; idx[0]<sz_x; ++idx[0])
- for (idx[1]=0; idx[1]<sz_y; ++idx[1])
- for (idx[2]=0; idx[2]<sz_z; ++idx[2])
+ itk::Index<4> idx4;
+ for (idx4[0]=0; idx4[0]<sz_x; ++idx4[0])
+ for (idx4[1]=0; idx4[1]<sz_y; ++idx4[1])
+ for (idx4[2]=0; idx4[2]<sz_z; ++idx4[2])
{
+ itk::Index<3> idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2];
+ if (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0)
+ continue;
+
vnl_vector_fixed<float,3> peak_dir;
vnl_vector_fixed<float,3> fitted_dir;
vnl_vector_fixed<float,3> overshoot_dir;
- for (idx[3]=0; idx[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx[3])
+ for (idx4[3]=0; idx4[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx4[3])
{
- peak_dir[idx[3]%3] = m_PeakImage->GetPixel(idx);
- fitted_dir[idx[3]%3] = m_FittedImage->GetPixel(idx);
- m_ResidualImage->SetPixel(idx, m_PeakImage->GetPixel(idx) - m_FittedImage->GetPixel(idx));
+ peak_dir[idx4[3]%3] = m_PeakImage->GetPixel(idx4);
+ fitted_dir[idx4[3]%3] = m_FittedImage->GetPixel(idx4);
+ m_ResidualImage->SetPixel(idx4, m_PeakImage->GetPixel(idx4) - m_FittedImage->GetPixel(idx4));
- if (idx[3]%3==2)
+ if (idx4[3]%3==2)
{
FD += peak_dir.magnitude();
- itk::Index<4> tidx= idx;
+ itk::Index<4> tidx= idx4;
if (peak_dir.magnitude()>fitted_dir.magnitude())
{
m_Coverage += fitted_dir.magnitude();
m_UnderexplainedImage->SetPixel(tidx, peak_dir[2]-fitted_dir[2]); tidx[3]--;
m_UnderexplainedImage->SetPixel(tidx, peak_dir[1]-fitted_dir[1]); tidx[3]--;
m_UnderexplainedImage->SetPixel(tidx, peak_dir[0]-fitted_dir[0]);
}
else
{
overshoot_dir[0] = fitted_dir[0]-peak_dir[0];
overshoot_dir[1] = fitted_dir[1]-peak_dir[1];
overshoot_dir[2] = fitted_dir[2]-peak_dir[2];
m_Coverage += peak_dir.magnitude();
m_Overshoot += overshoot_dir.magnitude();
m_OverexplainedImage->SetPixel(tidx, overshoot_dir[2]); tidx[3]--;
m_OverexplainedImage->SetPixel(tidx, overshoot_dir[1]); tidx[3]--;
m_OverexplainedImage->SetPixel(tidx, overshoot_dir[0]);
}
}
}
}
m_Coverage = m_Coverage/FD;
m_Overshoot = m_Overshoot/FD;
MITK_INFO << std::fixed << "Coverage: " << setprecision(1) << 100.0*m_Coverage << "%";
MITK_INFO << std::fixed << "Overshoot: " << setprecision(1) << 100.0*m_Overshoot << "%";
}
vnl_vector_fixed<float,3> FitFibersToImageFilter::GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer peak_image , vnl_vector_fixed<float,3> fiber_dir, int& id, double& w )
{
int m_NumDirs = peak_image->GetLargestPossibleRegion().GetSize()[3]/3;
vnl_vector_fixed<float,3> out_dir; out_dir.fill(0);
float angle = 0.8;
for (int i=0; i<m_NumDirs; i++)
{
vnl_vector_fixed<float,3> dir;
idx[3] = i*3;
dir[0] = peak_image->GetPixel(idx);
idx[3] += 1;
dir[1] = peak_image->GetPixel(idx);
idx[3] += 1;
dir[2] = peak_image->GetPixel(idx);
float mag = dir.magnitude();
if (mag<mitk::eps)
continue;
dir.normalize();
float a = dot_product(dir, fiber_dir);
if (fabs(a)>angle)
{
angle = fabs(a);
w = angle;
if (a<0)
out_dir = -dir;
else
out_dir = dir;
out_dir *= mag;
id = i;
}
}
return out_dir;
}
std::vector<mitk::FiberBundle::Pointer> FitFibersToImageFilter::GetTractograms() const
{
return m_Tractograms;
}
void FitFibersToImageFilter::SetTractograms(const std::vector<mitk::FiberBundle::Pointer> &tractograms)
{
m_Tractograms = tractograms;
}
}
diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h
index d784f3b2b7..119d39e57d 100644
--- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h
+++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h
@@ -1,252 +1,262 @@
#ifndef __itkFitFibersToImageFilter_h__
#define __itkFitFibersToImageFilter_h__
// MITK
#include <mitkFiberBundle.h>
#include <itkImageSource.h>
#include <mitkPeakImage.h>
#include <vnl/algo/vnl_lbfgsb.h>
#include <vnl/vnl_sparse_matrix.h>
#include <vnl/vnl_sparse_matrix_linear_system.h>
#include <itkImageDuplicator.h>
#include <itkTimeProbe.h>
#include <itkMersenneTwisterRandomVariateGenerator.h>
namespace itk{
/**
* \brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092). */
class FitFibersToImageFilter : public ImageSource< mitk::PeakImage::ItkPeakImageType >
{
public:
typedef FitFibersToImageFilter Self;
typedef ProcessObject Superclass;
typedef SmartPointer< Self > Pointer;
typedef SmartPointer< const Self > ConstPointer;
typedef itk::Point<float, 4> PointType4;
typedef mitk::PeakImage::ItkPeakImageType PeakImgType;
+ typedef itk::Image<unsigned char, 3> UcharImgType;
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
itkTypeMacro( FitFibersToImageFilter, ImageSource )
itkSetMacro( PeakImage, PeakImgType::Pointer)
itkGetMacro( PeakImage, PeakImgType::Pointer)
+ itkSetMacro( MaskImage, UcharImgType::Pointer)
+ itkGetMacro( MaskImage, UcharImgType::Pointer)
itkSetMacro( FitIndividualFibers, bool)
itkGetMacro( FitIndividualFibers, bool)
itkSetMacro( GradientTolerance, double)
itkGetMacro( GradientTolerance, double)
itkSetMacro( Lambda, double)
itkGetMacro( Lambda, double)
itkSetMacro( MaxIterations, int)
itkGetMacro( MaxIterations, int)
itkSetMacro( FiberSampling, float)
itkGetMacro( FiberSampling, float)
itkSetMacro( FilterOutliers, bool)
itkGetMacro( FilterOutliers, bool)
itkSetMacro( Verbose, bool)
itkGetMacro( Verbose, bool)
itkSetMacro( DeepCopy, bool)
itkGetMacro( DeepCopy, bool)
itkGetMacro( Weights, vnl_vector<double>)
itkGetMacro( FittedImage, PeakImgType::Pointer)
itkGetMacro( ResidualImage, PeakImgType::Pointer)
itkGetMacro( OverexplainedImage, PeakImgType::Pointer)
itkGetMacro( UnderexplainedImage, PeakImgType::Pointer)
itkGetMacro( Coverage, double)
itkGetMacro( Overshoot, double)
itkGetMacro( MeanWeight, double)
itkGetMacro( MedianWeight, double)
itkGetMacro( MinWeight, double)
itkGetMacro( MaxWeight, double)
+ itkGetMacro( NumUnknowns, unsigned int)
+ itkGetMacro( NumResiduals, unsigned int)
+ itkGetMacro( NumCoveredDirections, unsigned int)
void SetTractograms(const std::vector<mitk::FiberBundle::Pointer> &tractograms);
void GenerateData() override;
std::vector<mitk::FiberBundle::Pointer> GetTractograms() const;
protected:
FitFibersToImageFilter();
virtual ~FitFibersToImageFilter();
vnl_vector_fixed<float,3> GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer m_PeakImage , vnl_vector_fixed<float,3> fiber_dir, int& id, double& w );
std::vector< mitk::FiberBundle::Pointer > m_Tractograms;
PeakImgType::Pointer m_PeakImage;
+ UcharImgType::Pointer m_MaskImage;
bool m_FitIndividualFibers;
double m_GradientTolerance;
double m_Lambda;
int m_MaxIterations;
float m_FiberSampling;
double m_Coverage;
double m_Overshoot;
bool m_FilterOutliers;
double m_MeanWeight;
double m_MedianWeight;
double m_MinWeight;
double m_MaxWeight;
bool m_Verbose;
bool m_DeepCopy;
+ unsigned int m_NumUnknowns;
+ unsigned int m_NumResiduals;
+ unsigned int m_NumCoveredDirections;
// output
vnl_vector<double> m_Weights;
PeakImgType::Pointer m_UnderexplainedImage;
PeakImgType::Pointer m_OverexplainedImage;
PeakImgType::Pointer m_ResidualImage;
PeakImgType::Pointer m_FittedImage;
};
}
class VnlCostFunction : public vnl_cost_function
{
public:
vnl_sparse_matrix_linear_system< double >* S;
vnl_sparse_matrix< double > m_A;
vnl_sparse_matrix< double > m_A_Ones; // matrix indicating active weights with 1
vnl_vector< double > m_b;
double m_Lambda; // regularization factor
vnl_vector<double> row_sums; // number of active weights per row
vnl_vector<double> local_weight_means; // mean weight of each row
void SetProblem(vnl_sparse_matrix< double >& A, vnl_vector<double>& b, double lambda)
{
S = new vnl_sparse_matrix_linear_system<double>(A, b);
m_A = A;
m_b = b;
m_Lambda = lambda;
m_A_Ones.set_size(m_A.rows(), m_A.cols());
m_A.reset();
while (m_A.next())
m_A_Ones.put(m_A.getrow(), m_A.getcolumn(), 1);
unsigned int N = m_b.size();
vnl_vector<double> ones; ones.set_size(dim); ones.fill(1.0);
row_sums.set_size(N);
m_A_Ones.mult(ones, row_sums);
local_weight_means.set_size(N);
}
VnlCostFunction(const int NumVars) : vnl_cost_function(NumVars)
{
}
void regu_MSE(vnl_vector<double> const &x, double& cost)
{
double mean = x.mean();
vnl_vector<double> tx = x-mean;
cost += m_Lambda*1e8*tx.squared_magnitude()/x.size();
}
void regu_MSM(vnl_vector<double> const &x, double& cost)
{
cost += m_Lambda*1e8*x.squared_magnitude()/x.size();
}
void regu_localMSE(vnl_vector<double> const &x, double& cost)
{
m_A_Ones.mult(x, local_weight_means);
local_weight_means = element_quotient(local_weight_means, row_sums);
m_A_Ones.reset();
double regu = 0;
while (m_A_Ones.next())
{
double d = 0;
if (x[m_A_Ones.getcolumn()]>local_weight_means[m_A_Ones.getrow()])
d = std::exp(x[m_A_Ones.getcolumn()]) - std::exp(local_weight_means[m_A_Ones.getrow()]);
else
d = x[m_A_Ones.getcolumn()] - local_weight_means[m_A_Ones.getrow()];
regu += d*d;
}
cost += m_Lambda*regu/dim;
}
void grad_regu_MSE(vnl_vector<double> const &x, vnl_vector<double> &dx)
{
double mean = x.mean();
vnl_vector<double> tx = x-mean;
vnl_vector<double> tx2(dim, 0.0);
vnl_vector<double> h(dim, 1.0);
for (int c=0; c<dim; c++)
{
h[c] = dim-1;
tx2[c] += dot_product(h,tx);
h[c] = 1;
}
dx += tx2*m_Lambda*1e8*2.0/(dim*dim);
}
void grad_regu_MSM(vnl_vector<double> const &x, vnl_vector<double> &dx)
{
dx += m_Lambda*1e8*2.0*x/dim;
}
void grad_regu_localMSE(vnl_vector<double> const &x, vnl_vector<double> &dx)
{
m_A_Ones.mult(x, local_weight_means);
local_weight_means = element_quotient(local_weight_means, row_sums);
vnl_vector<double> exp_x = x.apply(std::exp);
vnl_vector<double> exp_means = local_weight_means.apply(std::exp);
vnl_vector<double> tdx(dim, 0);
m_A_Ones.reset();
while (m_A_Ones.next())
{
int c = m_A_Ones.getcolumn();
int r = m_A_Ones.getrow();
if (x[c]>local_weight_means[r])
tdx[c] += exp_x[c] * ( exp_x[c] - exp_means[r] );
else
tdx[c] += x[c] - local_weight_means[r];
}
dx += tdx*2.0*m_Lambda/dim;
}
double f(vnl_vector<double> const &x)
{
double cost = S->get_rms_error(x);
cost *= cost;
regu_localMSE(x, cost);
return cost;
}
void gradf(vnl_vector<double> const &x, vnl_vector<double> &dx)
{
dx.fill(0.0);
unsigned int N = m_b.size();
vnl_vector<double> d; d.set_size(N);
S->multiply(x,d);
d -= m_b;
S->transpose_multiply(d, dx);
dx *= 2.0/N;
grad_regu_localMSE(x,dx);
}
};
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFitFibersToImageFilter.cpp"
#endif
#endif // __itkFitFibersToImageFilter_h__
diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
index d53e7b4d7e..d33bf2cf32 100755
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
@@ -1,2325 +1,2335 @@
/*===================================================================
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, bool do_resampling)
{
vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
+ mitk::FiberBundle::Pointer fibCopy = this;
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 = this->GetDeepCopy();
fibCopy->ResampleLinear(minSpacing/5);
PolyData = fibCopy->GetFiberPolyData();
}
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
+ std::vector< float > new_weights;
+
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 (container->GetNumberOfPoints()>0)
+ {
+ new_weights.push_back(fibCopy->GetFiberWeight(i));
+ 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 newfib = mitk::FiberBundle::New(newPolyData);
+ for (unsigned int i=0; i<new_weights.size(); i++)
+ newfib->SetFiberWeight(i, new_weights.at(i));
+ return newfib;
}
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->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->GetNumberOfValues(); i++)
m_FiberWeights->SetValue(i, newWeight);
}
void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer<vtkFloatArray> weights)
{
if (m_NumFibers!=weights->GetNumberOfValues())
{
MITK_INFO << "Weights array not equal to number of fibers! " << weights->GetNumberOfValues() << " vs " << m_NumFibers;
return;
}
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/cmdapps/TractographyEvaluation/TractPlausibility.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp
index 899d63b403..ecf520ae8f 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp
@@ -1,201 +1,208 @@
/*===================================================================
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>
+#include <fstream>
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 > get_file_list(const std::string& path, const std::string& skipped)
{
+ std::srand(std::time(0));
+
std::vector< MaskType > mask_list;
itk::Directory::Pointer dir = itk::Directory::New();
+ ofstream myfile;
+ myfile.open (skipped);
if (dir->Load(path.c_str()))
{
int n = dir->GetNumberOfFiles();
for (int r = 0; r < n; r++)
{
const char *filename = dir->GetFile(r);
if (std::rand()%3==0 && ist::GetFilenameWithoutExtension(filename)!="CC")
{
MITK_INFO << "Dismissing " << ist::GetFilenameWithoutExtension(filename);
+ myfile << ist::GetFilenameName(filename) << "\n";
continue;
}
std::string ext = ist::GetFilenameExtension(filename);
if (ext==".nii" || ext==".nii.gz" || ext==".nrrd")
{
MITK_INFO << "Loading " << ist::GetFilenameWithoutExtension(filename);
streambuf *old = cout.rdbuf(); // <-- save
stringstream ss;
std::cout.rdbuf (ss.rdbuf()); // <-- redirect
MaskType m(LoadItkMaskImage(path + '/' + filename), ist::GetFilenameWithoutExtension(filename));
mask_list.push_back(m);
std::cout.rdbuf (old); // <-- restore
}
}
}
+ myfile.close();
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("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output folder", 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, "", "", 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"]);
+ string out_folder = us::any_cast<string>(parsedArgs["out"]);
try
{
- CreateFolderStructure(outRoot);
+ CreateFolderStructure(out_folder);
- std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder);
+ std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder, out_folder + "skipped_masks.txt");
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");
+ mitk::IOUtil::Save(not_gm_fibers, out_folder + "/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");
+ mitk::FiberBundle::Pointer known_tract = inputTractogram->ExtractFiberSubset(mask_image, true, false, false, 0.7, false);
+ mitk::IOUtil::Save(known_tract, out_folder + "/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_mrtrixspace.fib");
+ mitk::IOUtil::Save(all_known_tracts, out_folder + "/known_tracts/all_known_tracts.trk");
+ mitk::IOUtil::Save(TransformToMRtrixSpace(all_known_tracts), out_folder + "/known_tracts/all_known_tracts_mrtrixspace.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_mrtrixspace.fib");
+ mitk::IOUtil::Save(remaining_tracts, out_folder + "/plausible_tracts/remaining_tracts.trk");
+ mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), out_folder + "/plausible_tracts/remaining_tracts_mrtrixspace.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/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
index bceac56af1..6b1f1240c0 100755
--- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp
@@ -1,277 +1,292 @@
/*===================================================================
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 <itksys/SystemTools.hxx>
#include <itkDirectory.h>
#include <mitkFiberBundle.h>
#include <mitkPreferenceListReaderOptionsFunctor.h>
#include <itkImageFileWriter.h>
#include <mitkPeakImage.h>
#include <itkFitFibersToImageFilter.h>
#include <boost/lexical_cast.hpp>
using namespace std;
typedef itksys::SystemTools ist;
typedef itk::Point<float, 4> PointType4;
typedef itk::Image< float, 4 > PeakImgType;
std::vector< string > get_file_list(const std::string& path)
{
std::vector< string > file_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==".fib" || ext==".trk")
file_list.push_back(path + '/' + filename);
}
}
return file_list;
}
/*!
\brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092).
*/
int main(int argc, char* argv[])
{
mitkCommandLineParser parser;
parser.setTitle("");
parser.setCategory("Fiber Tracking Evaluation");
parser.setDescription("");
parser.setContributor("MIC");
parser.setArgumentPrefix("--", "-");
parser.addArgument("", "i1", mitkCommandLineParser::InputFile, "Input tractogram:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false);
parser.addArgument("", "i2", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false);
parser.addArgument("", "i3", mitkCommandLineParser::InputFile, "", "", us::Any(), false);
+ parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any());
+
+ parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "", "", us::Any(), false);
parser.addArgument("min_gain", "", mitkCommandLineParser::Float, "Min. gain:", "process stops if remaining bundles don't contribute enough", 0.05);
- parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false);
parser.addArgument("max_iter", "", mitkCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20);
parser.addArgument("bundle_based", "", mitkCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false);
parser.addArgument("min_g", "", mitkCommandLineParser::Float, "Min. gradient:", "lower termination threshold for gradient magnitude", 1e-5);
parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1);
parser.addArgument("dont_filter_outliers", "", mitkCommandLineParser::Bool, "Don't filter outliers:", "don't perform second optimization run with an upper weight bound based on the first weight estimation (95% quantile)", false);
map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
if (parsedArgs.size()==0)
return EXIT_FAILURE;
string fib_file = us::any_cast<string>(parsedArgs["i1"]);
string peak_file_name = us::any_cast<string>(parsedArgs["i2"]);
string candidate_folder = us::any_cast<string>(parsedArgs["i3"]);
string outRoot = us::any_cast<string>(parsedArgs["o"]);
bool single_fib = true;
if (parsedArgs.count("bundle_based"))
single_fib = !us::any_cast<bool>(parsedArgs["bundle_based"]);
int max_iter = 20;
if (parsedArgs.count("max_iter"))
max_iter = us::any_cast<int>(parsedArgs["max_iter"]);
float g_tol = 1e-5;
if (parsedArgs.count("min_g"))
g_tol = us::any_cast<float>(parsedArgs["min_g"]);
float min_gain = 0.05;
if (parsedArgs.count("min_gain"))
min_gain = us::any_cast<float>(parsedArgs["min_gain"]);
float lambda = 0.1;
if (parsedArgs.count("lambda"))
lambda = us::any_cast<float>(parsedArgs["lambda"]);
bool filter_outliers = true;
if (parsedArgs.count("dont_filter_outliers"))
filter_outliers = !us::any_cast<bool>(parsedArgs["dont_filter_outliers"]);
+ string mask_file = "";
+ if (parsedArgs.count("mask"))
+ mask_file = us::any_cast<string>(parsedArgs["mask"]);
+
try
{
itk::TimeProbe clock;
clock.Start();
mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {});
mitk::Image::Pointer inputImage = dynamic_cast<mitk::PeakImage*>(mitk::IOUtil::Load(peak_file_name, &functor)[0].GetPointer());
+ itk::FitFibersToImageFilter::UcharImgType::Pointer mask = nullptr;
+ if (mask_file.compare("")!=0)
+ {
+ mitk::Image::Pointer mitk_mask = dynamic_cast<mitk::Image*>(mitk::IOUtil::Load(mask_file)[0].GetPointer());
+ mitk::CastToItkImage(mitk_mask, mask);
+ }
+
typedef mitk::ImageToItk< PeakImgType > CasterType;
CasterType::Pointer caster = CasterType::New();
caster->SetInput(inputImage);
caster->Update();
PeakImgType::Pointer peak_image = caster->GetOutput();
std::vector< mitk::FiberBundle::Pointer > input_reference;
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(fib_file)[0].GetPointer());
if (fib.IsNull())
return EXIT_FAILURE;
input_reference.push_back(fib);
std::vector< mitk::FiberBundle::Pointer > input_candidates;
std::vector< string > candidate_tract_files = get_file_list(candidate_folder);
for (string f : candidate_tract_files)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(mitk::IOUtil::Load(f)[0].GetPointer());
if (fib.IsNull())
continue;
input_candidates.push_back(fib);
}
int iteration = 0;
std::string name = ist::GetFilenameWithoutExtension(fib_file);
itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
fitter->SetTractograms(input_reference);
fitter->SetFitIndividualFibers(single_fib);
fitter->SetMaxIterations(max_iter);
fitter->SetGradientTolerance(g_tol);
fitter->SetLambda(lambda);
fitter->SetFilterOutliers(filter_outliers);
fitter->SetPeakImage(peak_image);
fitter->SetVerbose(true);
fitter->SetDeepCopy(false);
+ fitter->SetMaskImage(mask);
fitter->Update();
// fitter->GetTractograms().at(0)->SetFiberWeights(fitter->GetCoverage());
// fitter->GetTractograms().at(0)->ColorFibersByFiberWeights(false, false);
mitk::IOUtil::Save(fitter->GetTractograms().at(0), outRoot + "0_" + name + ".fib");
peak_image = fitter->GetUnderexplainedImage();
itk::ImageFileWriter< PeakImgType >::Pointer writer = itk::ImageFileWriter< PeakImgType >::New();
writer->SetInput(peak_image);
writer->SetFileName(outRoot + boost::lexical_cast<string>(iteration) + "_" + name + ".nrrd");
writer->Update();
double coverage = fitter->GetCoverage();
MITK_INFO << "Iteration: " << iteration;
- MITK_INFO << std::fixed << "Coverage: " << setprecision(1) << 100.0*coverage << "%";
+ MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%";
// fitter->SetPeakImage(peak_image);
while (!input_candidates.empty())
{
streambuf *old = cout.rdbuf(); // <-- save
stringstream ss;
std::cout.rdbuf (ss.rdbuf()); // <-- redirect
double next_coverage = 0;
mitk::FiberBundle::Pointer best_candidate = nullptr;
for (auto fib : input_candidates)
{
// WHY NECESSARY AGAIN??
itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New();
fitter->SetFitIndividualFibers(single_fib);
fitter->SetMaxIterations(max_iter);
fitter->SetGradientTolerance(g_tol);
fitter->SetLambda(lambda);
fitter->SetFilterOutliers(filter_outliers);
fitter->SetVerbose(false);
fitter->SetPeakImage(peak_image);
fitter->SetDeepCopy(false);
+ fitter->SetMaskImage(mask);
// ******************************
fitter->SetTractograms({fib});
fitter->Update();
double candidate_coverage = fitter->GetCoverage();
if (candidate_coverage>next_coverage)
{
next_coverage = candidate_coverage;
if ((1.0-coverage) * next_coverage >= min_gain)
{
best_candidate = fitter->GetTractograms().at(0);
peak_image = fitter->GetUnderexplainedImage();
}
}
}
if (best_candidate.IsNull())
{
std::cout.rdbuf (old); // <-- restore
break;
}
// fitter->SetPeakImage(peak_image);
// best_candidate->SetFiberWeights((1.0-coverage) * next_coverage);
// best_candidate->ColorFibersByFiberWeights(false, false);
coverage += (1.0-coverage) * next_coverage;
int i=0;
std::vector< mitk::FiberBundle::Pointer > remaining_candidates;
std::vector< string > remaining_candidate_files;
for (auto fib : input_candidates)
{
if (fib!=best_candidate)
{
remaining_candidates.push_back(fib);
remaining_candidate_files.push_back(candidate_tract_files.at(i));
}
else
name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i));
++i;
}
input_candidates = remaining_candidates;
candidate_tract_files = remaining_candidate_files;
iteration++;
mitk::IOUtil::Save(best_candidate, outRoot + boost::lexical_cast<string>(iteration) + "_" + name + ".fib");
writer->SetInput(peak_image);
writer->SetFileName(outRoot + boost::lexical_cast<string>(iteration) + "_" + name + ".nrrd");
writer->Update();
std::cout.rdbuf (old); // <-- restore
MITK_INFO << "Iteration: " << iteration;
- MITK_INFO << std::fixed << "Coverage: " << setprecision(1) << 100.0*coverage << "% (+" << 100*(1.0-coverage) * next_coverage << "%)";
+ MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "% (+" << 100*(1.0-coverage) * next_coverage << "%)";
}
clock.Stop();
int h = clock.GetTotal()/3600;
int m = ((int)clock.GetTotal()%3600)/60;
int s = (int)clock.GetTotal()%60;
MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s";
}
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;
}