diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
index 9418e0878b..3941659fe4 100755
--- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp
@@ -1,2328 +1,2328 @@
/*===================================================================
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->InsertTupleValue(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 opacity, 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]);
if (opacity)
rgba[3] = (unsigned char) (255.0 * dev);
else
rgba[3] = (unsigned char) (255.0);
m_FiberColors->InsertTupleValue(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->InsertTupleValue(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->InsertTupleValue(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->InsertTupleValue(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)
{
vtkSmartPointer<vtkPolyData> PolyData = m_FiberPolyData;
if (anyPoint)
{
float minSpacing = 1;
if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
minSpacing = mask->GetSpacing()[0];
else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
minSpacing = mask->GetSpacing()[1];
else
minSpacing = mask->GetSpacing()[2];
mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy();
fibCopy->ResampleLinear(minSpacing/5);
PolyData = fibCopy->GetFiberPolyData();
}
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
- MITK_INFO << "Extracting fibers";
+ MITK_INFO << "Extracting fibers with mask image";
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp;
vtkCell* cell = PolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkCell* cellOriginal = m_FiberPolyData->GetCell(i);
int numPointsOriginal = cellOriginal->GetNumberOfPoints();
vtkPoints* pointsOriginal = cellOriginal->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
if (numPoints>1 && numPointsOriginal)
{
if (anyPoint)
{
int inside = 0;
int outside = 0;
if (!invert)
{
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
itk::Point<float, 3> itkP;
itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
itk::Index<3> idx;
mask->TransformPhysicalPointToIndex(itkP, idx);
if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 )
{
inside++;
if (fraction==0)
break;
}
else
outside++;
}
float current_fraction = 0.0;
if (inside+outside>0)
current_fraction = (float)inside/(inside+outside);
if (current_fraction>fraction)
{
for (int k=0; k<numPoints; k++)
{
double* p = points->GetPoint(k);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
else
{
bool includeFiber = true;
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
itk::Point<float, 3> itkP;
itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
itk::Index<3> idx;
mask->TransformPhysicalPointToIndex(itkP, idx);
if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) )
{
inside++;
includeFiber = false;
break;
}
else
outside++;
}
if (includeFiber)
{
for (int k=0; k<numPoints; k++)
{
double* p = points->GetPoint(k);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
}
else
{
double* start = pointsOriginal->GetPoint(0);
itk::Point<float, 3> itkStart;
itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2];
itk::Index<3> idxStart;
mask->TransformPhysicalPointToIndex(itkStart, idxStart);
double* end = pointsOriginal->GetPoint(numPointsOriginal-1);
itk::Point<float, 3> itkEnd;
itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2];
itk::Index<3> idxEnd;
mask->TransformPhysicalPointToIndex(itkEnd, idxEnd);
if (invert)
{
if (bothEnds)
{
if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 )
{
for (int j=0; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 )
{
for (int j=0; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
else
{
if (bothEnds)
{
if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) )
{
for (int j=0; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) )
{
for (int j=0; j<numPointsOriginal; j++)
{
double* p = pointsOriginal->GetPoint(j);
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
}
}
}
}
}
vtkNewCells->InsertNextCell(container);
}
if (vtkNewCells->GetNumberOfCells()<=0)
return nullptr;
vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
newPolyData->SetPoints(vtkNewPoints);
newPolyData->SetLines(vtkNewCells);
return mitk::FiberBundle::New(newPolyData);
}
mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert)
{
float minSpacing = 1;
if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
minSpacing = mask->GetSpacing()[0];
else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
minSpacing = mask->GetSpacing()[1];
else
minSpacing = mask->GetSpacing()[2];
mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy();
fibCopy->ResampleLinear(minSpacing/10);
vtkSmartPointer<vtkPolyData> PolyData =fibCopy->GetFiberPolyData();
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
MITK_INFO << "Cutting fibers";
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp;
vtkCell* cell = PolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
if (numPoints>1)
{
int newNumPoints = 0;
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
itk::Point<float, 3> itkP;
itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
itk::Index<3> idx;
mask->TransformPhysicalPointToIndex(itkP, idx);
if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert )
{
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
newNumPoints++;
}
else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert )
{
vtkIdType id = vtkNewPoints->InsertNextPoint(p);
container->GetPointIds()->InsertNextId(id);
newNumPoints++;
}
else if (newNumPoints>0)
{
vtkNewCells->InsertNextCell(container);
newNumPoints = 0;
container = vtkSmartPointer<vtkPolyLine>::New();
}
}
if (newNumPoints>0)
vtkNewCells->InsertNextCell(container);
}
}
if (vtkNewCells->GetNumberOfCells()<=0)
return nullptr;
vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
newPolyData->SetPoints(vtkNewPoints);
newPolyData->SetLines(vtkNewCells);
mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData);
newFib->Compress(0.1);
return newFib;
}
mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage)
{
if (roi==nullptr || !(dynamic_cast<PlanarFigure*>(roi->GetData()) || dynamic_cast<PlanarFigureComposite*>(roi->GetData())) )
return nullptr;
std::vector<long> tmp = ExtractFiberIdSubset(roi, storage);
if (tmp.size()<=0)
return mitk::FiberBundle::New();
vtkSmartPointer<vtkPolyData> pTmp = GeneratePolyDataByIds(tmp);
return mitk::FiberBundle::New(pTmp);
}
std::vector<long> mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage)
{
std::vector<long> result;
if (roi==nullptr || roi->GetData()==nullptr)
return result;
mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast<mitk::PlanarFigureComposite*>(roi->GetData());
if (!pfc.IsNull()) // handle composite
{
DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi);
if (children->size()==0)
return result;
switch (pfc->getOperationType())
{
case 0: // AND
{
MITK_INFO << "AND";
result = this->ExtractFiberIdSubset(children->ElementAt(0), storage);
std::vector<long>::iterator it;
for (unsigned int i=1; i<children->Size(); ++i)
{
std::vector<long> inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage);
std::vector<long> rest(std::min(result.size(),inRoi.size()));
it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
rest.resize( it - rest.begin() );
result = rest;
}
break;
}
case 1: // OR
{
MITK_INFO << "OR";
result = ExtractFiberIdSubset(children->ElementAt(0), storage);
std::vector<long>::iterator it;
for (unsigned int i=1; i<children->Size(); ++i)
{
it = result.end();
std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
result.insert(it, inRoi.begin(), inRoi.end());
}
// remove duplicates
sort(result.begin(), result.end());
it = unique(result.begin(), result.end());
result.resize( it - result.begin() );
break;
}
case 2: // NOT
{
MITK_INFO << "NOT";
for(long i=0; i<this->GetNumFibers(); i++)
result.push_back(i);
std::vector<long>::iterator it;
for (unsigned int i=0; i<children->Size(); ++i)
{
std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
std::vector<long> rest(result.size()-inRoi.size());
it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
rest.resize( it - rest.begin() );
result = rest;
}
break;
}
}
}
else if ( dynamic_cast<mitk::PlanarFigure*>(roi->GetData()) ) // actual extraction
{
if ( dynamic_cast<mitk::PlanarPolygon*>(roi->GetData()) )
{
mitk::PlanarFigure::Pointer planarPoly = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
//create vtkPolygon using controlpoints from planarFigure polygon
vtkSmartPointer<vtkPolygon> polygonVtk = vtkSmartPointer<vtkPolygon>::New();
for (unsigned int i=0; i<planarPoly->GetNumberOfControlPoints(); ++i)
{
itk::Point<double,3> p = planarPoly->GetWorldControlPoint(i);
vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] );
polygonVtk->GetPointIds()->InsertNextId(id);
}
MITK_INFO << "Extracting with polygon";
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp ;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
for (int j=0; j<numPoints-1; j++)
{
// Inputs
double p1[3] = {0,0,0};
points->GetPoint(j, p1);
double p2[3] = {0,0,0};
points->GetPoint(j+1, p2);
double tolerance = 0.001;
// Outputs
double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
double x[3] = {0,0,0}; // The coordinate of the intersection
double pcoords[3] = {0,0,0};
int subId = 0;
int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId);
if (iD!=0)
{
result.push_back(i);
break;
}
}
}
}
else if ( dynamic_cast<mitk::PlanarCircle*>(roi->GetData()) )
{
mitk::PlanarFigure::Pointer planarFigure = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal();
planeNormal.Normalize();
//calculate circle radius
mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint
mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint
double radius = V1w.EuclideanDistanceTo(V2w);
radius *= radius;
MITK_INFO << "Extracting with circle";
boost::progress_display disp(m_NumFibers);
for (int i=0; i<m_NumFibers; i++)
{
++disp ;
vtkCell* cell = m_FiberPolyData->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
for (int j=0; j<numPoints-1; j++)
{
// Inputs
double p1[3] = {0,0,0};
points->GetPoint(j, p1);
double p2[3] = {0,0,0};
points->GetPoint(j+1, p2);
// Outputs
double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
double x[3] = {0,0,0}; // The coordinate of the intersection
int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x);
if (iD!=0)
{
double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]);
if( dist <= radius)
{
result.push_back(i);
break;
}
}
}
}
}
return result;
}
return result;
}
void mitk::FiberBundle::UpdateFiberGeometry()
{
vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
cleaner->SetInputData(m_FiberPolyData);
cleaner->PointMergingOff();
cleaner->Update();
m_FiberPolyData = cleaner->GetOutput();
m_FiberLengths.clear();
m_MeanFiberLength = 0;
m_MedianFiberLength = 0;
m_LengthStDev = 0;
m_NumFibers = m_FiberPolyData->GetNumberOfCells();
if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints())
this->ColorFibersByOrientation();
if (m_FiberWeights->GetSize()!=m_NumFibers)
{
m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
m_FiberWeights->SetName("FIBER_WEIGHTS");
m_FiberWeights->SetNumberOfValues(m_NumFibers);
this->SetFiberWeights(1);
}
if (m_NumFibers<=0) // no fibers present; apply default geometry
{
m_MinFiberLength = 0;
m_MaxFiberLength = 0;
mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
geometry->SetImageGeometry(false);
float b[] = {0, 1, 0, 1, 0, 1};
geometry->SetFloatBounds(b);
SetGeometry(geometry);
return;
}
double b[6];
m_FiberPolyData->GetBounds(b);
// calculate statistics
for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
{
vtkCell* cell = m_FiberPolyData->GetCell(i);
int p = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
float length = 0;
for (int j=0; j<p-1; j++)
{
double p1[3];
points->GetPoint(j, p1);
double p2[3];
points->GetPoint(j+1, p2);
float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
length += dist;
}
m_FiberLengths.push_back(length);
m_MeanFiberLength += length;
if (i==0)
{
m_MinFiberLength = length;
m_MaxFiberLength = length;
}
else
{
if (length<m_MinFiberLength)
m_MinFiberLength = length;
if (length>m_MaxFiberLength)
m_MaxFiberLength = length;
}
}
m_MeanFiberLength /= m_NumFibers;
std::vector< float > sortedLengths = m_FiberLengths;
std::sort(sortedLengths.begin(), sortedLengths.end());
for (int i=0; i<m_NumFibers; i++)
m_LengthStDev += (m_MeanFiberLength-sortedLengths.at(i))*(m_MeanFiberLength-sortedLengths.at(i));
if (m_NumFibers>1)
m_LengthStDev /= (m_NumFibers-1);
else
m_LengthStDev = 0;
m_LengthStDev = std::sqrt(m_LengthStDev);
m_MedianFiberLength = sortedLengths.at(m_NumFibers/2);
mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
geometry->SetFloatBounds(b);
this->SetGeometry(geometry);
m_UpdateTime3D.Modified();
m_UpdateTime2D.Modified();
}
float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const
{
return m_FiberWeights->GetValue(fiber);
}
void mitk::FiberBundle::SetFiberWeights(float newWeight)
{
for (int i=0; i<m_FiberWeights->GetSize(); i++)
m_FiberWeights->SetValue(i, newWeight);
}
void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer<vtkFloatArray> weights)
{
if (m_NumFibers!=weights->GetSize())
{
MITK_INFO << "Weights array not equal to number of fibers!";
return;
}
for (int i=0; i<weights->GetSize(); 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->GetTupleValue(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->InsertTupleValue(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/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp
index 3b85ee12d7..5c897e61b1 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp
@@ -1,1533 +1,1540 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
// Blueberry
#include <berryISelectionService.h>
#include <berryIWorkbenchPart.h>
#include <berryIWorkbenchWindow.h>
// Qmitk
#include "QmitkFiberProcessingView.h"
// Qt
#include <QMessageBox>
// MITK
#include <mitkNodePredicateProperty.h>
#include <mitkImageCast.h>
#include <mitkPointSet.h>
#include <mitkPlanarCircle.h>
#include <mitkPlanarPolygon.h>
#include <mitkPlanarRectangle.h>
#include <mitkPlanarFigureInteractor.h>
#include <mitkImageAccessByItk.h>
#include <mitkDataNodeObject.h>
#include <mitkTensorImage.h>
#include "usModuleRegistry.h"
#include <itkFiberCurvatureFilter.h>
#include "mitkNodePredicateDataType.h"
#include <mitkNodePredicateProperty.h>
#include <mitkNodePredicateAnd.h>
#include <mitkNodePredicateNot.h>
#include <mitkNodePredicateOr.h>
// ITK
#include <itkResampleImageFilter.h>
#include <itkGaussianInterpolateImageFunction.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkTractsToFiberEndingsImageFilter.h>
#include <itkTractDensityImageFilter.h>
#include <itkImageRegion.h>
#include <itkTractsToRgbaImageFilter.h>
#define _USE_MATH_DEFINES
#include <math.h>
const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing";
const std::string id_DataManager = "org.mitk.views.datamanager";
using namespace mitk;
QmitkFiberProcessingView::QmitkFiberProcessingView()
: QmitkAbstractView()
, m_Controls( 0 )
, m_CircleCounter(0)
, m_PolygonCounter(0)
, m_UpsamplingFactor(1)
{
}
// Destructor
QmitkFiberProcessingView::~QmitkFiberProcessingView()
{
RemoveObservers();
}
void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent )
{
// build up qt view, unless already done
if ( !m_Controls )
{
// create GUI widgets from the Qt Designer's .ui file
m_Controls = new Ui::QmitkFiberProcessingViewControls;
m_Controls->setupUi( parent );
connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) );
connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) );
connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) );
connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) );
connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) );
connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) );
connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) );
connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) );
connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) );
connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract()));
connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove()));
connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify()));
connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged()));
m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage());
mitk::TNodePredicateDataType<mitk::Image>::Pointer isMitkImage = mitk::TNodePredicateDataType<mitk::Image>::New();
mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage");
mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage");
mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage");
mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti);
isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi);
mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage);
mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage);
m_Controls->m_ColorMapBox->SetPredicate(finalPredicate);
m_Controls->label_17->setVisible(false);
m_Controls->m_FiberExtractionFractionBox->setVisible(false);
}
UpdateGui();
}
void QmitkFiberProcessingView::OnMaskExtractionChanged()
{
if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2)
{
m_Controls->label_17->setVisible(true);
m_Controls->m_FiberExtractionFractionBox->setVisible(true);
m_Controls->m_BothEnds->setVisible(false);
}
else
{
m_Controls->label_17->setVisible(false);
m_Controls->m_FiberExtractionFractionBox->setVisible(false);
if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3)
m_Controls->m_BothEnds->setVisible(true);
}
}
void QmitkFiberProcessingView::SetFocus()
{
m_Controls->toolBoxx->setFocus();
}
void QmitkFiberProcessingView::Modify()
{
switch (m_Controls->m_ModificationMethodBox->currentIndex())
{
case 0:
{
ResampleSelectedBundlesSpline();
break;
}
case 1:
{
ResampleSelectedBundlesLinear();
break;
}
case 2:
{
CompressSelectedBundles();
break;
}
case 3:
{
DoImageColorCoding();
break;
}
case 4:
{
MirrorFibers();
break;
}
case 5:
{
WeightFibers();
break;
}
case 6:
{
DoCurvatureColorCoding();
break;
}
case 7:
{
DoWeightColorCoding();
break;
}
}
}
void QmitkFiberProcessingView::WeightFibers()
{
float weight = this->m_Controls->m_BundleWeightBox->value();
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
fib->SetFiberWeights(weight);
}
}
void QmitkFiberProcessingView::Remove()
{
switch (m_Controls->m_RemovalMethodBox->currentIndex())
{
case 0:
{
RemoveDir();
break;
}
case 1:
{
PruneBundle();
break;
}
case 2:
{
ApplyCurvatureThreshold();
break;
}
case 3:
{
RemoveWithMask(false);
break;
}
case 4:
{
RemoveWithMask(true);
break;
}
}
}
void QmitkFiberProcessingView::Extract()
{
switch (m_Controls->m_ExtractionMethodBox->currentIndex())
{
case 0:
{
ExtractWithPlanarFigure();
break;
}
case 1:
{
switch (m_Controls->m_ExtractionBoxMask->currentIndex())
{
{
case 0:
ExtractWithMask(true, false);
break;
}
{
case 1:
ExtractWithMask(true, true);
break;
}
{
case 2:
ExtractWithMask(false, false);
break;
}
{
case 3:
ExtractWithMask(false, true);
break;
}
}
break;
}
}
}
void QmitkFiberProcessingView::PruneBundle()
{
int minLength = this->m_Controls->m_PruneFibersMinBox->value();
int maxLength = this->m_Controls->m_PruneFibersMaxBox->value();
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
if (!fib->RemoveShortFibers(minLength))
QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
else if (!fib->RemoveLongFibers(maxLength))
QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
}
RenderingManager::GetInstance()->RequestUpdateAll();
}
void QmitkFiberProcessingView::ApplyCurvatureThreshold()
{
int angle = this->m_Controls->m_CurvSpinBox->value();
int dist = this->m_Controls->m_CurvDistanceSpinBox->value();
std::vector< DataNode::Pointer > nodes = m_SelectedFB;
for (auto node : nodes)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New();
filter->SetInputFiberBundle(fib);
filter->SetAngularDeviation(angle);
filter->SetDistance(dist);
filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked());
filter->Update();
mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle();
if (newFib->GetNumFibers()>0)
{
newFib->ColorFibersByOrientation();
node->SetData(newFib);
}
else
QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
}
RenderingManager::GetInstance()->RequestUpdateAll();
}
void QmitkFiberProcessingView::RemoveDir()
{
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
vnl_vector_fixed<double,3> dir;
dir[0] = m_Controls->m_ExtractDirX->value();
dir[1] = m_Controls->m_ExtractDirY->value();
dir[2] = m_Controls->m_ExtractDirZ->value();
fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*M_PI/180));
}
RenderingManager::GetInstance()->RequestUpdateAll();
}
void QmitkFiberProcessingView::RemoveWithMask(bool removeInside)
{
if (m_MaskImageNode.IsNull())
return;
mitk::Image::Pointer mitkMask = dynamic_cast<mitk::Image*>(m_MaskImageNode->GetData());
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
itkUCharImageType::Pointer mask = itkUCharImageType::New();
mitk::CastToItkImage(mitkMask, mask);
mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside);
if (newFib->GetNumFibers()<=0)
{
QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
continue;
}
node->SetData(newFib);
}
RenderingManager::GetInstance()->RequestUpdateAll();
}
void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert)
{
if (m_MaskImageNode.IsNull())
return;
mitk::Image::Pointer mitkMask = dynamic_cast<mitk::Image*>(m_MaskImageNode->GetData());
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
QString name(node->GetName().c_str());
itkUCharImageType::Pointer mask = itkUCharImageType::New();
mitk::CastToItkImage(mitkMask, mask);
mitk::FiberBundle::Pointer newFib = fib->ExtractFiberSubset(mask, !onlyEnds, invert, m_Controls->m_BothEnds->isChecked(), m_Controls->m_FiberExtractionFractionBox->value());
if (newFib->GetNumFibers()<=0)
{
QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
continue;
}
DataNode::Pointer newNode = DataNode::New();
newNode->SetData(newFib);
if (invert)
{
name += "_not";
if (onlyEnds)
name += "-ending-in-mask";
else
name += "-passing-mask";
}
else
{
if (onlyEnds)
name += "_ending-in-mask";
else
name += "_passing-mask";
}
newNode->SetName(name.toStdString());
GetDataStorage()->Add(newNode);
node->SetVisibility(false);
}
}
void QmitkFiberProcessingView::GenerateRoiImage()
{
if (m_SelectedPF.empty())
return;
mitk::BaseGeometry::Pointer geometry;
if (!m_SelectedFB.empty())
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.front()->GetData());
geometry = fib->GetGeometry();
}
else if (m_SelectedImage)
geometry = m_SelectedImage->GetGeometry();
else
return;
itk::Vector<double,3> spacing = geometry->GetSpacing();
spacing /= m_UpsamplingFactor;
mitk::Point3D newOrigin = geometry->GetOrigin();
mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds();
newOrigin[0] += bounds.GetElement(0);
newOrigin[1] += bounds.GetElement(2);
newOrigin[2] += bounds.GetElement(4);
itk::Matrix<double, 3, 3> direction;
itk::ImageRegion<3> imageRegion;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j];
imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor);
imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor);
imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor);
m_PlanarFigureImage = itkUCharImageType::New();
m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing
m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin
m_PlanarFigureImage->SetDirection( direction ); // Set the image direction
m_PlanarFigureImage->SetRegions( imageRegion );
m_PlanarFigureImage->Allocate();
m_PlanarFigureImage->FillBuffer( 0 );
Image::Pointer tmpImage = Image::New();
tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer());
tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer());
std::string name = m_SelectedPF.at(0)->GetName();
WritePfToImage(m_SelectedPF.at(0), tmpImage);
for (unsigned int i=1; i<m_SelectedPF.size(); i++)
{
name += "+";
name += m_SelectedPF.at(i)->GetName();
WritePfToImage(m_SelectedPF.at(i), tmpImage);
}
DataNode::Pointer node = DataNode::New();
tmpImage = Image::New();
tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer());
tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer());
node->SetData(tmpImage);
node->SetName(name);
this->GetDataStorage()->Add(node);
}
void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image)
{
if (dynamic_cast<mitk::PlanarFigure*>(node->GetData()))
{
m_PlanarFigure = dynamic_cast<mitk::PlanarFigure*>(node->GetData());
AccessFixedDimensionByItk_2(
image,
InternalReorientImagePlane, 3,
m_PlanarFigure->GetGeometry(), -1);
AccessFixedDimensionByItk_2(
m_InternalImage,
InternalCalculateMaskFromPlanarFigure,
3, 2, node->GetName() );
}
else if (dynamic_cast<mitk::PlanarFigureComposite*>(node->GetData()))
{
DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node);
for (unsigned int i=0; i<children->Size(); i++)
{
WritePfToImage(children->at(i), image);
}
}
}
template < typename TPixel, unsigned int VImageDimension >
void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex )
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::Image< float, VImageDimension > FloatImageType;
typedef itk::ResampleImageFilter<ImageType, FloatImageType, double> ResamplerType;
typename ResamplerType::Pointer resampler = ResamplerType::New();
mitk::PlaneGeometry* planegeo = dynamic_cast<mitk::PlaneGeometry*>(planegeo3D);
float upsamp = m_UpsamplingFactor;
float gausssigma = 0.5;
// Spacing
typename ResamplerType::SpacingType spacing = planegeo->GetSpacing();
spacing[0] = image->GetSpacing()[0] / upsamp;
spacing[1] = image->GetSpacing()[1] / upsamp;
spacing[2] = image->GetSpacing()[2];
resampler->SetOutputSpacing( spacing );
// Size
typename ResamplerType::SizeType size;
size[0] = planegeo->GetExtentInMM(0) / spacing[0];
size[1] = planegeo->GetExtentInMM(1) / spacing[1];
size[2] = 1;
resampler->SetSize( size );
// Origin
typename mitk::Point3D orig = planegeo->GetOrigin();
typename mitk::Point3D corrorig;
planegeo3D->WorldToIndex(orig,corrorig);
corrorig[0] += 0.5/upsamp;
corrorig[1] += 0.5/upsamp;
corrorig[2] += 0;
planegeo3D->IndexToWorld(corrorig,corrorig);
resampler->SetOutputOrigin(corrorig );
// Direction
typename ResamplerType::DirectionType direction;
typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix();
for(unsigned int c=0; c<matrix.ColumnDimensions; c++)
{
double sum = 0;
for(unsigned int r=0; r<matrix.RowDimensions; r++)
sum += matrix(r,c)*matrix(r,c);
for(unsigned int r=0; r<matrix.RowDimensions; r++)
direction(r,c) = matrix(r,c)/sqrt(sum);
}
resampler->SetOutputDirection( direction );
// Gaussian interpolation
if(gausssigma != 0)
{
double sigma[3];
for( unsigned int d = 0; d < 3; d++ )
sigma[d] = gausssigma * image->GetSpacing()[d];
double alpha = 2.0;
typedef itk::GaussianInterpolateImageFunction<ImageType, double> GaussianInterpolatorType;
typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New();
interpolator->SetInputImage( image );
interpolator->SetParameters( sigma, alpha );
resampler->SetInterpolator( interpolator );
}
else
{
typedef typename itk::LinearInterpolateImageFunction<ImageType, double> InterpolatorType;
typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
interpolator->SetInputImage( image );
resampler->SetInterpolator( interpolator );
}
resampler->SetInput( image );
resampler->SetDefaultPixelValue(0);
resampler->Update();
if(additionalIndex < 0)
{
this->m_InternalImage = mitk::Image::New();
this->m_InternalImage->InitializeByItk( resampler->GetOutput() );
this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() );
}
}
template < typename TPixel, unsigned int VImageDimension >
void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string )
{
typedef itk::Image< TPixel, VImageDimension > ImageType;
typedef itk::CastImageFilter< ImageType, itkUCharImageType > CastFilterType;
// Generate mask image as new image with same header as input image and
// initialize with "1".
itkUCharImageType::Pointer newMaskImage = itkUCharImageType::New();
newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing
newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin
newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction
newMaskImage->SetRegions( image->GetLargestPossibleRegion() );
newMaskImage->Allocate();
newMaskImage->FillBuffer( 1 );
// Generate VTK polygon from (closed) PlanarFigure polyline
// (The polyline points are shifted by -0.5 in z-direction to make sure
// that the extrusion filter, which afterwards elevates all points by +0.5
// in z-direction, creates a 3D object which is cut by the the plane z=0)
const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 );
vtkPolyData *polyline = vtkPolyData::New();
polyline->Allocate( 1, 1 );
// Determine x- and y-dimensions depending on principal axis
int i0, i1;
switch ( axis )
{
case 0:
i0 = 1;
i1 = 2;
break;
case 1:
i0 = 0;
i1 = 2;
break;
case 2:
default:
i0 = 0;
i1 = 1;
break;
}
// Create VTK polydata object of polyline contour
vtkPoints *points = vtkPoints::New();
PlanarFigure::PolyLineType::const_iterator it;
unsigned int numberOfPoints = 0;
for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it )
{
Point3D point3D;
// Convert 2D point back to the local index coordinates of the selected image
Point2D point2D = *it;
planarFigurePlaneGeometry->WorldToIndex(point2D, point2D);
point2D[0] -= 0.5/m_UpsamplingFactor;
point2D[1] -= 0.5/m_UpsamplingFactor;
planarFigurePlaneGeometry->IndexToWorld(point2D, point2D);
planarFigurePlaneGeometry->Map( point2D, point3D );
// Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil
if ( !imageGeometry3D->IsInside( point3D ) )
{
float bounds[2] = {0,0};
bounds[0] =
this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0);
bounds[1] =
this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1);
imageGeometry3D->WorldToIndex( point3D, point3D );
if (point3D[i0]<0)
point3D[i0] = 0.0;
else if (point3D[i0]>bounds[0])
point3D[i0] = bounds[0]-0.001;
if (point3D[i1]<0)
point3D[i1] = 0.0;
else if (point3D[i1]>bounds[1])
point3D[i1] = bounds[1]-0.001;
points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 );
numberOfPoints++;
}
else
{
imageGeometry3D->WorldToIndex( point3D, point3D );
// Add point to polyline array
points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 );
numberOfPoints++;
}
}
polyline->SetPoints( points );
points->Delete();
vtkIdType *ptIds = new vtkIdType[numberOfPoints];
for ( vtkIdType i = 0; i < numberOfPoints; ++i )
ptIds[i] = i;
polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds );
// Extrude the generated contour polygon
vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New();
extrudeFilter->SetInputData( polyline );
extrudeFilter->SetScaleFactor( 1 );
extrudeFilter->SetExtrusionTypeToNormalExtrusion();
extrudeFilter->SetVector( 0.0, 0.0, 1.0 );
// Make a stencil from the extruded polygon
vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New();
polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() );
// Export from ITK to VTK (to use a VTK filter)
typedef itk::VTKImageImport< itkUCharImageType > ImageImportType;
typedef itk::VTKImageExport< itkUCharImageType > ImageExportType;
typename ImageExportType::Pointer itkExporter = ImageExportType::New();
itkExporter->SetInput( newMaskImage );
vtkImageImport *vtkImporter = vtkImageImport::New();
this->ConnectPipelines( itkExporter, vtkImporter );
vtkImporter->Update();
// Apply the generated image stencil to the input image
vtkImageStencil *imageStencilFilter = vtkImageStencil::New();
imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() );
imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() );
imageStencilFilter->ReverseStencilOff();
imageStencilFilter->SetBackgroundValue( 0 );
imageStencilFilter->Update();
// Export from VTK back to ITK
vtkImageExport *vtkExporter = vtkImageExport::New();
vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() );
vtkExporter->Update();
typename ImageImportType::Pointer itkImporter = ImageImportType::New();
this->ConnectPipelines( vtkExporter, itkImporter );
itkImporter->Update();
// calculate cropping bounding box
m_InternalImageMask3D = itkImporter->GetOutput();
m_InternalImageMask3D->SetDirection(image->GetDirection());
itk::ImageRegionConstIterator<itkUCharImageType>
itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion());
itk::ImageRegionIterator<ImageType>
itimage(image, image->GetLargestPossibleRegion());
itmask.GoToBegin();
itimage.GoToBegin();
typename ImageType::SizeType lowersize = {{itk::NumericTraits<unsigned long>::max(),itk::NumericTraits<unsigned long>::max(),itk::NumericTraits<unsigned long>::max()}};
typename ImageType::SizeType uppersize = {{0,0,0}};
while( !itmask.IsAtEnd() )
{
if(itmask.Get() == 0)
itimage.Set(0);
else
{
typename ImageType::IndexType index = itimage.GetIndex();
typename ImageType::SizeType signedindex;
signedindex[0] = index[0];
signedindex[1] = index[1];
signedindex[2] = index[2];
lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0];
lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1];
lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2];
uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0];
uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1];
uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2];
}
++itmask;
++itimage;
}
typename ImageType::IndexType index;
index[0] = lowersize[0];
index[1] = lowersize[1];
index[2] = lowersize[2];
typename ImageType::SizeType size;
size[0] = uppersize[0] - lowersize[0] + 1;
size[1] = uppersize[1] - lowersize[1] + 1;
size[2] = uppersize[2] - lowersize[2] + 1;
itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size);
// crop internal mask
typedef itk::RegionOfInterestImageFilter< itkUCharImageType, itkUCharImageType > ROIMaskFilterType;
typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New();
roi2->SetRegionOfInterest(cropRegion);
roi2->SetInput(m_InternalImageMask3D);
roi2->Update();
m_InternalImageMask3D = roi2->GetOutput();
Image::Pointer tmpImage = Image::New();
tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer());
tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer());
Image::Pointer tmpImage2 = Image::New();
tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer());
const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 );
const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 );
typedef itk::ImageRegionIteratorWithIndex<itkUCharImageType> IteratorType;
IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion());
imageIterator.GoToBegin();
while ( !imageIterator.IsAtEnd() )
{
unsigned char val = imageIterator.Value();
if (val>0)
{
itk::Index<3> index = imageIterator.GetIndex();
Point3D point;
point[0] = index[0];
point[1] = index[1];
point[2] = index[2];
intImageGeometry3D->IndexToWorld(point, point);
pfImageGeometry3D->WorldToIndex(point, point);
point[i0] += 0.5;
point[i1] += 0.5;
index[0] = point[0];
index[1] = point[1];
index[2] = point[2];
if (pfImageGeometry3D->IsIndexInside(index))
m_PlanarFigureImage->SetPixel(index, 1);
}
++imageIterator;
}
// Clean up VTK objects
polyline->Delete();
extrudeFilter->Delete();
polyDataToImageStencil->Delete();
vtkImporter->Delete();
imageStencilFilter->Delete();
//vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak??
delete[] ptIds;
}
void QmitkFiberProcessingView::UpdateGui()
{
m_Controls->m_FibLabel->setText("<font color='red'>mandatory</font>");
m_Controls->m_PfLabel->setText("<font color='grey'>needed for extraction</font>");
m_Controls->m_InputData->setTitle("Please Select Input Data");
m_Controls->m_RemoveButton->setEnabled(false);
m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false);
m_Controls->PFCompoANDButton->setEnabled(false);
m_Controls->PFCompoORButton->setEnabled(false);
m_Controls->PFCompoNOTButton->setEnabled(false);
m_Controls->m_GenerateRoiImage->setEnabled(false);
m_Controls->m_ExtractFibersButton->setEnabled(false);
m_Controls->m_ModifyButton->setEnabled(false);
m_Controls->m_CopyBundle->setEnabled(false);
m_Controls->m_JoinBundles->setEnabled(false);
m_Controls->m_SubstractBundles->setEnabled(false);
// disable alle frames
m_Controls->m_BundleWeightFrame->setVisible(false);
m_Controls->m_ExtactionFramePF->setVisible(false);
m_Controls->m_RemoveDirectionFrame->setVisible(false);
m_Controls->m_RemoveLengthFrame->setVisible(false);
m_Controls->m_RemoveCurvatureFrame->setVisible(false);
m_Controls->m_SmoothFibersFrame->setVisible(false);
m_Controls->m_CompressFibersFrame->setVisible(false);
m_Controls->m_ColorFibersFrame->setVisible(false);
m_Controls->m_MirrorFibersFrame->setVisible(false);
m_Controls->m_MaskExtractionFrame->setVisible(false);
m_Controls->m_ColorMapBox->setVisible(false);
bool pfSelected = !m_SelectedPF.empty();
bool fibSelected = !m_SelectedFB.empty();
bool multipleFibsSelected = (m_SelectedFB.size()>1);
bool maskSelected = m_MaskImageNode.IsNotNull();
bool imageSelected = m_SelectedImage.IsNotNull();
// toggle visibility of elements according to selected method
switch ( m_Controls->m_ExtractionMethodBox->currentIndex() )
{
case 0:
m_Controls->m_ExtactionFramePF->setVisible(true);
break;
case 1:
m_Controls->m_MaskExtractionFrame->setVisible(true);
break;
}
switch ( m_Controls->m_RemovalMethodBox->currentIndex() )
{
case 0:
m_Controls->m_RemoveDirectionFrame->setVisible(true);
if ( fibSelected )
m_Controls->m_RemoveButton->setEnabled(true);
break;
case 1:
m_Controls->m_RemoveLengthFrame->setVisible(true);
if ( fibSelected )
m_Controls->m_RemoveButton->setEnabled(true);
break;
case 2:
m_Controls->m_RemoveCurvatureFrame->setVisible(true);
if ( fibSelected )
m_Controls->m_RemoveButton->setEnabled(true);
break;
case 3:
break;
case 4:
break;
}
switch ( m_Controls->m_ModificationMethodBox->currentIndex() )
{
case 0:
m_Controls->m_SmoothFibersFrame->setVisible(true);
break;
case 1:
m_Controls->m_SmoothFibersFrame->setVisible(true);
break;
case 2:
m_Controls->m_CompressFibersFrame->setVisible(true);
break;
case 3:
m_Controls->m_ColorFibersFrame->setVisible(true);
m_Controls->m_ColorMapBox->setVisible(true);
break;
case 4:
m_Controls->m_MirrorFibersFrame->setVisible(true);
if (m_SelectedSurfaces.size()>0)
m_Controls->m_ModifyButton->setEnabled(true);
break;
case 5:
m_Controls->m_BundleWeightFrame->setVisible(true);
break;
case 6:
m_Controls->m_ColorFibersFrame->setVisible(true);
break;
case 7:
m_Controls->m_ColorFibersFrame->setVisible(true);
break;
}
// are fiber bundles selected?
if ( fibSelected )
{
m_Controls->m_CopyBundle->setEnabled(true);
m_Controls->m_ModifyButton->setEnabled(true);
m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true);
m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str()));
// one bundle and one planar figure needed to extract fibers
if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0)
{
m_Controls->m_InputData->setTitle("Input Data");
m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str()));
m_Controls->m_ExtractFibersButton->setEnabled(true);
}
// more than two bundles needed to join/subtract
if (multipleFibsSelected)
{
m_Controls->m_FibLabel->setText("multiple bundles selected");
m_Controls->m_JoinBundles->setEnabled(true);
m_Controls->m_SubstractBundles->setEnabled(true);
}
if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1)
{
m_Controls->m_InputData->setTitle("Input Data");
m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str()));
m_Controls->m_ExtractFibersButton->setEnabled(true);
}
if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) )
{
m_Controls->m_InputData->setTitle("Input Data");
m_Controls->m_PfLabel->setText(QString(m_MaskImageNode->GetName().c_str()));
m_Controls->m_RemoveButton->setEnabled(true);
}
}
// are planar figures selected?
if (pfSelected)
{
if ( fibSelected || m_SelectedImage.IsNotNull())
m_Controls->m_GenerateRoiImage->setEnabled(true);
if (m_SelectedPF.size() > 1)
{
m_Controls->PFCompoANDButton->setEnabled(true);
m_Controls->PFCompoORButton->setEnabled(true);
}
else
m_Controls->PFCompoNOTButton->setEnabled(true);
}
// is image selected
if (imageSelected || maskSelected)
{
m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true);
}
}
-void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* )
+void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node )
{
+ for (auto fnode: m_SelectedFB)
+ if (node == fnode)
+ {
+ m_SelectedFB.clear();
+ break;
+ }
+
berry::IWorkbenchPart::Pointer nullPart;
QList<mitk::DataNode::Pointer> nodes;
OnSelectionChanged(nullPart, nodes);
}
void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* )
{
if (!m_Controls->m_InteractiveBox->isChecked())
{
berry::IWorkbenchPart::Pointer nullPart;
QList<mitk::DataNode::Pointer> nodes;
OnSelectionChanged(nullPart, nodes);
}
}
-void QmitkFiberProcessingView::OnStartInteraction()
-{
-
-}
-
void QmitkFiberProcessingView::OnEndInteraction()
{
if (m_Controls->m_InteractiveBox->isChecked())
ExtractWithPlanarFigure(true);
}
void QmitkFiberProcessingView::AddObservers()
{
typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType;
for (auto node : m_SelectedPF)
{
mitk::PlanarFigure* figure = dynamic_cast<mitk::PlanarFigure*>(node->GetData());
if (figure!=nullptr)
{
- // add observer for event when interaction with figure starts
- SimpleCommandType::Pointer startInteractionCommand = SimpleCommandType::New();
- startInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnStartInteraction);
- m_StartInteractionObserverTag = figure->AddObserver( mitk::StartInteractionPlanarFigureEvent(), startInteractionCommand );
+ figure->RemoveAllObservers();
// add observer for event when interaction with figure starts
SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New();
endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction);
m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand );
}
}
}
void QmitkFiberProcessingView::RemoveObservers()
{
for (auto node : m_SelectedPF)
{
mitk::PlanarFigure* figure = dynamic_cast<mitk::PlanarFigure*>(node->GetData());
if (figure!=nullptr)
- {
- // figure->RemoveAllObservers();
- if (figure->HasObserver(mitk::StartInteractionPlanarFigureEvent()))
- figure->RemoveObserver(m_StartInteractionObserverTag);
- if (figure->HasObserver(mitk::EndInteractionPlanarFigureEvent()))
- figure->RemoveObserver(m_EndInteractionObserverTag);
- m_StartInteractionObserverTag = -1;
- m_EndInteractionObserverTag = -1;
- }
+ figure->RemoveAllObservers();
}
}
void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList<mitk::DataNode::Pointer>& nodes)
{
RemoveObservers();
//reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection
+ std::vector<mitk::DataNode::Pointer> lastSelectedFB = m_SelectedFB;
m_SelectedFB.clear();
m_SelectedPF.clear();
m_SelectedSurfaces.clear();
m_SelectedImage = nullptr;
m_MaskImageNode = nullptr;
for (auto node: nodes)
{
if ( dynamic_cast<mitk::FiberBundle*>(node->GetData()) )
m_SelectedFB.push_back(node);
else if (dynamic_cast<mitk::PlanarPolygon*>(node->GetData()) || dynamic_cast<mitk::PlanarFigureComposite*>(node->GetData()) || dynamic_cast<mitk::PlanarCircle*>(node->GetData()))
m_SelectedPF.push_back(node);
else if (dynamic_cast<mitk::Image*>(node->GetData()))
{
m_SelectedImage = dynamic_cast<mitk::Image*>(node->GetData());
bool isBinary = false;
node->GetPropertyValue<bool>("binary", isBinary);
if (isBinary)
m_MaskImageNode = node;
}
else if (dynamic_cast<mitk::Surface*>(node->GetData()))
m_SelectedSurfaces.push_back(dynamic_cast<mitk::Surface*>(node->GetData()));
}
+ // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle
+ if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked())
+ m_SelectedFB = lastSelectedFB;
+
+ // if no fibers or surfaces are selected, select topmost
if (m_SelectedFB.empty() && m_SelectedSurfaces.empty())
{
int maxLayer = 0;
itk::VectorContainer<unsigned int, mitk::DataNode::Pointer>::ConstPointer nodes = this->GetDataStorage()->GetAll();
for (unsigned int i=0; i<nodes->Size(); i++)
if (dynamic_cast<mitk::FiberBundle*>(nodes->at(i)->GetData()))
{
mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i));
if (sources->Size()>0)
continue;
int layer = 0;
nodes->at(i)->GetPropertyValue("layer", layer);
if (layer>=maxLayer)
{
maxLayer = layer;
m_SelectedFB.clear();
m_SelectedFB.push_back(nodes->at(i));
}
}
}
+ // if no plar figure is selected, select topmost
if (m_SelectedPF.empty())
{
int maxLayer = 0;
itk::VectorContainer<unsigned int, mitk::DataNode::Pointer>::ConstPointer nodes = this->GetDataStorage()->GetAll();
for (unsigned int i=0; i<nodes->Size(); i++)
if (dynamic_cast<mitk::PlanarPolygon*>(nodes->at(i)->GetData()) || dynamic_cast<mitk::PlanarFigureComposite*>(nodes->at(i)->GetData()) || dynamic_cast<mitk::PlanarCircle*>(nodes->at(i)->GetData()))
{
mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i));
if (sources->Size()>0)
continue;
int layer = 0;
nodes->at(i)->GetPropertyValue("layer", layer);
if (layer>=maxLayer)
{
maxLayer = layer;
m_SelectedPF.clear();
m_SelectedPF.push_back(nodes->at(i));
}
}
}
AddObservers();
UpdateGui();
}
void QmitkFiberProcessingView::OnDrawPolygon()
{
mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New();
figure->ClosedOn();
this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter));
}
void QmitkFiberProcessingView::OnDrawCircle()
{
mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New();
this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter));
}
void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* )
{
// initialize figure's geometry with empty geometry
mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New();
figure->SetPlaneGeometry( emptygeometry );
//set desired data to DataNode where Planarfigure is stored
mitk::DataNode::Pointer newNode = mitk::DataNode::New();
newNode->SetName(name.toStdString());
newNode->SetData(figure);
newNode->SetBoolProperty("planarfigure.3drendering", true);
newNode->SetBoolProperty("planarfigure.3drendering.fill", true);
mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast<mitk::PlanarFigureInteractor*>(newNode->GetDataInteractor().GetPointer());
if(figureInteractor.IsNull())
{
figureInteractor = mitk::PlanarFigureInteractor::New();
us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" );
figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule );
figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule );
figureInteractor->SetDataNode(newNode);
}
// figure drawn on the topmost layer / image
GetDataStorage()->Add(newNode );
+ RemoveObservers();
for(unsigned int i = 0; i < m_SelectedPF.size(); i++)
m_SelectedPF[i]->SetSelected(false);
newNode->SetSelected(true);
m_SelectedPF.clear();
m_SelectedPF.push_back(newNode);
+ AddObservers();
UpdateGui();
}
void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive)
{
if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){
QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!");
return;
}
std::vector<mitk::DataNode::Pointer> fiberBundles = m_SelectedFB;
mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0);
for (unsigned int i=0; i<fiberBundles.size(); i++)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(fiberBundles.at(i)->GetData());
mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage());
if (interactive && m_Controls->m_InteractiveBox->isChecked())
{
if (extFB->GetNumFibers()<=0)
continue;
if (m_InteractiveNode.IsNull())
{
m_InteractiveNode = mitk::DataNode::New();
QString name("Interactive");
m_InteractiveNode->SetName(name.toStdString());
GetDataStorage()->Add(m_InteractiveNode);
}
+ float op = 5.0/sqrt(fib->GetNumFibers());
+
fib->SetFiberColors(255, 255, 255);
- fiberBundles.at(i)->SetFloatProperty("opacity", 0.01);
+ fiberBundles.at(i)->SetFloatProperty("opacity", op);
fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false);
m_InteractiveNode->SetData(extFB);
if (auto renderWindowPart = this->GetRenderWindowPart())
renderWindowPart->RequestUpdate();
}
else
{
if (extFB->GetNumFibers()<=0)
{
QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers.");
continue;
}
mitk::DataNode::Pointer node;
node = mitk::DataNode::New();
node->SetData(extFB);
QString name(fiberBundles.at(i)->GetName().c_str());
name += "*";
node->SetName(name.toStdString());
fiberBundles.at(i)->SetVisibility(false);
GetDataStorage()->Add(node);
}
}
}
void QmitkFiberProcessingView::GenerateAndComposite()
{
mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New();
PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND);
mitk::DataNode::Pointer newPFCNode;
newPFCNode = mitk::DataNode::New();
newPFCNode->SetName("AND");
newPFCNode->SetData(PFCAnd);
AddCompositeToDatastorage(newPFCNode, m_SelectedPF);
+ RemoveObservers();
m_SelectedPF.clear();
m_SelectedPF.push_back(newPFCNode);
+ AddObservers();
UpdateGui();
}
void QmitkFiberProcessingView::GenerateOrComposite()
{
mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New();
PFCOr->setOperationType(mitk::PlanarFigureComposite::OR);
mitk::DataNode::Pointer newPFCNode;
newPFCNode = mitk::DataNode::New();
newPFCNode->SetName("OR");
newPFCNode->SetData(PFCOr);
+ RemoveObservers();
AddCompositeToDatastorage(newPFCNode, m_SelectedPF);
m_SelectedPF.clear();
m_SelectedPF.push_back(newPFCNode);
UpdateGui();
}
void QmitkFiberProcessingView::GenerateNotComposite()
{
mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New();
PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT);
mitk::DataNode::Pointer newPFCNode;
newPFCNode = mitk::DataNode::New();
newPFCNode->SetName("NOT");
newPFCNode->SetData(PFCNot);
+ RemoveObservers();
AddCompositeToDatastorage(newPFCNode, m_SelectedPF);
m_SelectedPF.clear();
m_SelectedPF.push_back(newPFCNode);
+ AddObservers();
UpdateGui();
}
void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector<mitk::DataNode::Pointer> children, mitk::DataNode::Pointer parentNode )
{
pfc->SetSelected(true);
if (parentNode.IsNotNull())
GetDataStorage()->Add(pfc, parentNode);
else
GetDataStorage()->Add(pfc);
for (auto child : children)
{
if (dynamic_cast<PlanarFigure*>(child->GetData()))
{
mitk::DataNode::Pointer newChild;
newChild = mitk::DataNode::New();
newChild->SetData(dynamic_cast<PlanarFigure*>(child->GetData()));
newChild->SetName( child->GetName() );
newChild->SetBoolProperty("planarfigure.3drendering", true);
newChild->SetBoolProperty("planarfigure.3drendering.fill", true);
GetDataStorage()->Add(newChild, pfc);
GetDataStorage()->Remove(child);
}
else if (dynamic_cast<PlanarFigureComposite*>(child->GetData()))
{
mitk::DataNode::Pointer newChild;
newChild = mitk::DataNode::New();
newChild->SetData(dynamic_cast<PlanarFigureComposite*>(child->GetData()));
newChild->SetName( child->GetName() );
std::vector< mitk::DataNode::Pointer > grandChildVector;
mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child);
for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it )
grandChildVector.push_back(*it);
AddCompositeToDatastorage(newChild, grandChildVector, pfc);
GetDataStorage()->Remove(child);
}
}
UpdateGui();
}
void QmitkFiberProcessingView::CopyBundles()
{
if ( m_SelectedFB.empty() ){
QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!");
MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!";
return;
}
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy();
node->SetVisibility(false);
QString name("");
name += QString(m_SelectedFB.at(0)->GetName().c_str());
name += "_copy";
mitk::DataNode::Pointer fbNode = mitk::DataNode::New();
fbNode->SetData(newFib);
fbNode->SetName(name.toStdString());
fbNode->SetVisibility(true);
GetDataStorage()->Add(fbNode);
}
UpdateGui();
}
void QmitkFiberProcessingView::JoinBundles()
{
if ( m_SelectedFB.size()<2 ){
QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!");
MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!";
return;
}
mitk::FiberBundle::Pointer newBundle = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(0)->GetData());
m_SelectedFB.at(0)->SetVisibility(false);
QString name("");
name += QString(m_SelectedFB.at(0)->GetName().c_str());
for (unsigned int i=1; i<m_SelectedFB.size(); i++)
{
newBundle = newBundle->AddBundle(dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(i)->GetData()));
name += "+"+QString(m_SelectedFB.at(i)->GetName().c_str());
m_SelectedFB.at(i)->SetVisibility(false);
}
mitk::DataNode::Pointer fbNode = mitk::DataNode::New();
fbNode->SetData(newBundle);
fbNode->SetName(name.toStdString());
fbNode->SetVisibility(true);
GetDataStorage()->Add(fbNode);
UpdateGui();
}
void QmitkFiberProcessingView::SubstractBundles()
{
if ( m_SelectedFB.size()<2 ){
QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!");
MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!";
return;
}
mitk::FiberBundle::Pointer newBundle = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(0)->GetData());
m_SelectedFB.at(0)->SetVisibility(false);
QString name("");
name += QString(m_SelectedFB.at(0)->GetName().c_str());
for (unsigned int i=1; i<m_SelectedFB.size(); i++)
{
newBundle = newBundle->SubtractBundle(dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.at(i)->GetData()));
if (newBundle.IsNull())
break;
name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str());
m_SelectedFB.at(i)->SetVisibility(false);
}
if (newBundle.IsNull())
{
QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!");
return;
}
mitk::DataNode::Pointer fbNode = mitk::DataNode::New();
fbNode->SetData(newBundle);
fbNode->SetName(name.toStdString());
fbNode->SetVisibility(true);
GetDataStorage()->Add(fbNode);
UpdateGui();
}
void QmitkFiberProcessingView::ResampleSelectedBundlesSpline()
{
double factor = this->m_Controls->m_SmoothFibersBox->value();
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
fib->ResampleSpline(factor);
}
RenderingManager::GetInstance()->RequestUpdateAll();
}
void QmitkFiberProcessingView::ResampleSelectedBundlesLinear()
{
double factor = this->m_Controls->m_SmoothFibersBox->value();
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
fib->ResampleLinear(factor);
}
RenderingManager::GetInstance()->RequestUpdateAll();
}
void QmitkFiberProcessingView::CompressSelectedBundles()
{
double factor = this->m_Controls->m_ErrorThresholdBox->value();
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
fib->Compress(factor);
fib->ColorFibersByOrientation();
}
RenderingManager::GetInstance()->RequestUpdateAll();
}
void QmitkFiberProcessingView::DoImageColorCoding()
{
if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull())
{
QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available.");
return;
}
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
fib->ColorFibersByScalarMap(dynamic_cast<mitk::Image*>(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked());
}
if (auto renderWindowPart = this->GetRenderWindowPart())
{
renderWindowPart->RequestUpdate();
}
}
void QmitkFiberProcessingView::DoCurvatureColorCoding()
{
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked());
}
if (auto renderWindowPart = this->GetRenderWindowPart())
{
renderWindowPart->RequestUpdate();
}
}
void QmitkFiberProcessingView::DoWeightColorCoding()
{
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked());
}
if (auto renderWindowPart = this->GetRenderWindowPart())
{
renderWindowPart->RequestUpdate();
}
}
void QmitkFiberProcessingView::MirrorFibers()
{
unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex();
for (auto node : m_SelectedFB)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
if (m_SelectedImage.IsNotNull())
fib->SetReferenceGeometry(m_SelectedImage->GetGeometry());
fib->MirrorFibers(axis);
}
for (auto surf : m_SelectedSurfaces)
{
vtkSmartPointer<vtkPolyData> poly = surf->GetVtkPolyData();
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
for (int i=0; i<poly->GetNumberOfPoints(); i++)
{
double* point = poly->GetPoint(i);
point[axis] *= -1;
vtkNewPoints->InsertNextPoint(point);
}
poly->SetPoints(vtkNewPoints);
surf->CalculateBoundingBox();
}
if (auto renderWindowPart = this->GetRenderWindowPart())
{
renderWindowPart->RequestUpdate();
}
}
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h
index b9aaaaf09e..522d304afe 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h
@@ -1,200 +1,199 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#ifndef QmitkFiberProcessingView_h
#define QmitkFiberProcessingView_h
#include <QmitkAbstractView.h>
#include "ui_QmitkFiberProcessingViewControls.h"
#include <mitkPlanarFigureComposite.h>
#include <mitkFiberBundle.h>
#include <mitkSurface.h>
#include <itkCastImageFilter.h>
#include <itkVTKImageImport.h>
#include <itkVTKImageExport.h>
#include <itkRegionOfInterestImageFilter.h>
#include <vtkLinearExtrusionFilter.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkSelectEnclosedPoints.h>
#include <vtkImageImport.h>
#include <vtkImageExport.h>
#include <vtkImageStencil.h>
#include <vtkSmartPointer.h>
#include <vtkSelection.h>
#include <vtkSelectionNode.h>
#include <vtkExtractSelectedThresholds.h>
#include <vtkFloatArray.h>
/*!
\brief View to process fiber bundles. Supplies methods to extract fibers from the bundle, fiber resampling, mirroring, join and subtract bundles and much more.
*/
class QmitkFiberProcessingView : public QmitkAbstractView
{
// this is needed for all Qt objects that should have a Qt meta-object
// (everything that derives from QObject and wants to have signal/slots)
Q_OBJECT
public:
typedef itk::Image< unsigned char, 3 > itkUCharImageType;
static const std::string VIEW_ID;
QmitkFiberProcessingView();
virtual ~QmitkFiberProcessingView();
virtual void CreateQtPartControl(QWidget *parent) override;
///
/// Sets the focus to an internal widget.
///
virtual void SetFocus() override;
protected slots:
void OnDrawCircle(); ///< add circle interactors etc.
void OnDrawPolygon(); ///< add circle interactors etc.
void GenerateAndComposite();
void GenerateOrComposite();
void GenerateNotComposite();
void CopyBundles(); ///< add copies of selected bundles to data storage
void JoinBundles(); ///< merge selected fiber bundles
void SubstractBundles(); ///< subtract bundle A from bundle B. Not commutative! Defined by order of selection.
void GenerateRoiImage(); ///< generate binary image of selected planar figures.
void Remove();
void Extract();
void Modify();
void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection
void OnMaskExtractionChanged();
virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey = nullptr, mitk::BaseProperty *property = nullptr );
protected:
void MirrorFibers(); ///< mirror bundle on the specified plane
void ResampleSelectedBundlesSpline(); ///<
void ResampleSelectedBundlesLinear(); ///<
void DoImageColorCoding(); ///< color fibers by selected scalar image
void DoWeightColorCoding(); ///< color fibers by their respective weights
void DoCurvatureColorCoding(); ///< color fibers by curvature
void CompressSelectedBundles(); ///< remove points below certain error threshold
void WeightFibers();
void RemoveWithMask(bool removeInside);
void RemoveDir();
void ApplyCurvatureThreshold(); ///< remove/split fibers with a too high curvature threshold
void PruneBundle(); ///< remove too short/too long fibers
void ExtractWithMask(bool onlyEnds, bool invert);
void ExtractWithPlanarFigure(bool interactive=false);
void OnEndInteraction();
- void OnStartInteraction();
/// \brief called by QmitkAbstractView when DataManager's selection has changed
virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList<mitk::DataNode::Pointer>& nodes) override;
Ui::QmitkFiberProcessingViewControls* m_Controls;
/** Connection from VTK to ITK */
template <typename VTK_Exporter, typename ITK_Importer>
void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer)
{
importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback());
importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback());
importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback());
importer->SetSpacingCallback(exporter->GetSpacingCallback());
importer->SetOriginCallback(exporter->GetOriginCallback());
importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback());
importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback());
importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback());
importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback());
importer->SetDataExtentCallback(exporter->GetDataExtentCallback());
importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback());
importer->SetCallbackUserData(exporter->GetCallbackUserData());
}
template <typename ITK_Exporter, typename VTK_Importer>
void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer)
{
importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback());
importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback());
importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback());
importer->SetSpacingCallback(exporter->GetSpacingCallback());
importer->SetOriginCallback(exporter->GetOriginCallback());
importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback());
importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback());
importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback());
importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback());
importer->SetDataExtentCallback(exporter->GetDataExtentCallback());
importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback());
importer->SetCallbackUserData(exporter->GetCallbackUserData());
}
template < typename TPixel, unsigned int VImageDimension >
void InternalCalculateMaskFromPlanarFigure(
itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string nodeName );
template < typename TPixel, unsigned int VImageDimension >
void InternalReorientImagePlane(
const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex );
int m_CircleCounter; ///< used for data node naming
int m_PolygonCounter; ///< used for data node naming
std::vector<mitk::DataNode::Pointer> m_SelectedFB; ///< selected fiber bundle nodes
std::vector<mitk::DataNode::Pointer> m_SelectedPF; ///< selected planar figure nodes
std::vector<mitk::Surface::Pointer> m_SelectedSurfaces;
mitk::Image::Pointer m_SelectedImage;
mitk::Image::Pointer m_InternalImage;
mitk::PlanarFigure::Pointer m_PlanarFigure;
itkUCharImageType::Pointer m_InternalImageMask3D;
itkUCharImageType::Pointer m_PlanarFigureImage;
float m_UpsamplingFactor; ///< upsampling factor for all image generations
mitk::DataNode::Pointer m_MaskImageNode;
unsigned int m_StartInteractionObserverTag;
unsigned int m_EndInteractionObserverTag;
mitk::DataNode::Pointer m_InteractiveNode;
void AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector<mitk::DataNode::Pointer> children, mitk::DataNode::Pointer parentNode=nullptr);
void debugPFComposition(mitk::PlanarFigureComposite::Pointer , int );
void WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image);
mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute);
mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundle::Pointer fib);
mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib);
mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib);
void NodeAdded( const mitk::DataNode* node ) override;
void NodeRemoved(const mitk::DataNode* node) override;
void RemoveObservers();
void AddObservers();
};
#endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED