Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
mitkFiberBundle.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 
18 #define _USE_MATH_DEFINES
19 #include "mitkFiberBundle.h"
20 
21 #include <mitkPlanarCircle.h>
22 #include <mitkPlanarPolygon.h>
23 #include <mitkPlanarFigureComposite.h>
24 #include "mitkImagePixelReadAccessor.h"
25 #include <mitkPixelTypeMultiplex.h>
26 
27 #include <vtkPointData.h>
28 #include <vtkDataArray.h>
29 #include <vtkUnsignedCharArray.h>
30 #include <vtkPolyLine.h>
31 #include <vtkCellArray.h>
32 #include <vtkCellData.h>
33 #include <vtkIdFilter.h>
34 #include <vtkClipPolyData.h>
35 #include <vtkPlane.h>
36 #include <vtkDoubleArray.h>
37 #include <vtkKochanekSpline.h>
38 #include <vtkParametricFunctionSource.h>
39 #include <vtkParametricSpline.h>
40 #include <vtkPolygon.h>
41 #include <vtkCleanPolyData.h>
42 #include <cmath>
43 #include <boost/progress.hpp>
44 #include <vtkTransformPolyDataFilter.h>
45 #include <mitkTransferFunction.h>
46 #include <vtkLookupTable.h>
47 #include <mitkLookupTable.h>
48 
49 const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs";
50 
51 using namespace std;
52 
53 mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData )
54  : m_NumFibers(0)
55  , m_FiberSampling(0)
56 {
57  m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
58  m_FiberWeights->SetName("FIBER_WEIGHTS");
59 
60  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
61  if (fiberPolyData != nullptr)
62  {
63  m_FiberPolyData = fiberPolyData;
64  this->ColorFibersByOrientation();
65  }
66 
67  this->UpdateFiberGeometry();
68  this->GenerateFiberIds();
69 }
70 
71 mitk::FiberBundle::~FiberBundle()
72 {
73 
74 }
75 
76 mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy()
77 {
78  mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData);
79  newFib->SetFiberColors(this->m_FiberColors);
80  newFib->SetFiberWeights(this->m_FiberWeights);
81  return newFib;
82 }
83 
84 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GeneratePolyDataByIds(std::vector<long> fiberIds)
85 {
86  vtkSmartPointer<vtkPolyData> newFiberPolyData = vtkSmartPointer<vtkPolyData>::New();
87  vtkSmartPointer<vtkCellArray> newLineSet = vtkSmartPointer<vtkCellArray>::New();
88  vtkSmartPointer<vtkPoints> newPointSet = vtkSmartPointer<vtkPoints>::New();
89 
90  auto finIt = fiberIds.begin();
91  while ( finIt != fiberIds.end() )
92  {
93  if (*finIt < 0 || *finIt>GetNumFibers()){
94  MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt;
95  break;
96  }
97 
98  vtkSmartPointer<vtkCell> fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber);
99  vtkSmartPointer<vtkPoints> fibPoints = fiber->GetPoints();
100  vtkSmartPointer<vtkPolyLine> newFiber = vtkSmartPointer<vtkPolyLine>::New();
101  newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() );
102 
103  for(int i=0; i<fibPoints->GetNumberOfPoints(); i++)
104  {
105  newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints());
106  newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]);
107  }
108 
109  newLineSet->InsertNextCell(newFiber);
110  ++finIt;
111  }
112 
113  newFiberPolyData->SetPoints(newPointSet);
114  newFiberPolyData->SetLines(newLineSet);
115  return newFiberPolyData;
116 }
117 
118 // merge two fiber bundles
119 mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib)
120 {
121  if (fib==nullptr)
122  {
123  MITK_WARN << "trying to call AddBundle with NULL argument";
124  return nullptr;
125  }
126  MITK_INFO << "Adding fibers";
127 
128  vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
129  vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
130  vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
131 
132  // add current fiber bundle
133  vtkSmartPointer<vtkFloatArray> weights = vtkSmartPointer<vtkFloatArray>::New();
134  weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers());
135 
136  unsigned int counter = 0;
137  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
138  {
139  vtkCell* cell = m_FiberPolyData->GetCell(i);
140  int numPoints = cell->GetNumberOfPoints();
141  vtkPoints* points = cell->GetPoints();
142 
143  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
144  for (int j=0; j<numPoints; j++)
145  {
146  double p[3];
147  points->GetPoint(j, p);
148 
149  vtkIdType id = vNewPoints->InsertNextPoint(p);
150  container->GetPointIds()->InsertNextId(id);
151  }
152  weights->InsertValue(counter, this->GetFiberWeight(i));
153  vNewLines->InsertNextCell(container);
154  counter++;
155  }
156 
157  // add new fiber bundle
158  for (int i=0; i<fib->GetFiberPolyData()->GetNumberOfCells(); i++)
159  {
160  vtkCell* cell = fib->GetFiberPolyData()->GetCell(i);
161  int numPoints = cell->GetNumberOfPoints();
162  vtkPoints* points = cell->GetPoints();
163 
164  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
165  for (int j=0; j<numPoints; j++)
166  {
167  double p[3];
168  points->GetPoint(j, p);
169 
170  vtkIdType id = vNewPoints->InsertNextPoint(p);
171  container->GetPointIds()->InsertNextId(id);
172  }
173  weights->InsertValue(counter, fib->GetFiberWeight(i));
174  vNewLines->InsertNextCell(container);
175  counter++;
176  }
177 
178  // initialize polydata
179  vNewPolyData->SetPoints(vNewPoints);
180  vNewPolyData->SetLines(vNewLines);
181 
182  // initialize fiber bundle
183  mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData);
184  newFib->SetFiberWeights(weights);
185  return newFib;
186 }
187 
188 // subtract two fiber bundles
189 mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib)
190 {
191  MITK_INFO << "Subtracting fibers";
192 
193  vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
194  vtkSmartPointer<vtkCellArray> vNewLines = vtkSmartPointer<vtkCellArray>::New();
195  vtkSmartPointer<vtkPoints> vNewPoints = vtkSmartPointer<vtkPoints>::New();
196 
197  // iterate over current fibers
198  boost::progress_display disp(m_NumFibers);
199  for( int i=0; i<m_NumFibers; i++ )
200  {
201  ++disp;
202  vtkCell* cell = m_FiberPolyData->GetCell(i);
203  int numPoints = cell->GetNumberOfPoints();
204  vtkPoints* points = cell->GetPoints();
205 
206  if (points==nullptr || numPoints<=0)
207  continue;
208 
209  int numFibers2 = fib->GetNumFibers();
210  bool contained = false;
211  for( int i2=0; i2<numFibers2; i2++ )
212  {
213  vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i2);
214  int numPoints2 = cell2->GetNumberOfPoints();
215  vtkPoints* points2 = cell2->GetPoints();
216 
217  if (points2==nullptr)// || numPoints2<=0)
218  continue;
219 
220  // check endpoints
221  if (numPoints2==numPoints)
222  {
223  itk::Point<float, 3> point_start = GetItkPoint(points->GetPoint(0));
224  itk::Point<float, 3> point_end = GetItkPoint(points->GetPoint(numPoints-1));
225  itk::Point<float, 3> point2_start = GetItkPoint(points2->GetPoint(0));
226  itk::Point<float, 3> point2_end = GetItkPoint(points2->GetPoint(numPoints2-1));
227 
228  if ((point_start.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps) ||
229  (point_start.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps))
230  {
231  // further checking ???
232  contained = true;
233  break;
234  }
235  }
236  }
237 
238  // add to result because fiber is not subtracted
239  if (!contained)
240  {
241  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
242  for( int j=0; j<numPoints; j++)
243  {
244  vtkIdType id = vNewPoints->InsertNextPoint(points->GetPoint(j));
245  container->GetPointIds()->InsertNextId(id);
246  }
247  vNewLines->InsertNextCell(container);
248  }
249  }
250  if(vNewLines->GetNumberOfCells()==0)
251  return nullptr;
252  // initialize polydata
253  vNewPolyData->SetPoints(vNewPoints);
254  vNewPolyData->SetLines(vNewLines);
255 
256  // initialize fiber bundle
257  return mitk::FiberBundle::New(vNewPolyData);
258 }
259 
260 itk::Point<float, 3> mitk::FiberBundle::GetItkPoint(double point[3])
261 {
262  itk::Point<float, 3> itkPoint;
263  itkPoint[0] = point[0];
264  itkPoint[1] = point[1];
265  itkPoint[2] = point[2];
266  return itkPoint;
267 }
268 
269 /*
270  * set polydata (additional flag to recompute fiber geometry, default = true)
271  */
272 void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer<vtkPolyData> fiberPD, bool updateGeometry)
273 {
274  if (fiberPD == nullptr)
275  this->m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
276  else
277  {
278  m_FiberPolyData->DeepCopy(fiberPD);
279  ColorFibersByOrientation();
280  }
281 
282  m_NumFibers = m_FiberPolyData->GetNumberOfLines();
283 
284  if (updateGeometry)
285  UpdateFiberGeometry();
286  GenerateFiberIds();
287 }
288 
289 /*
290  * return vtkPolyData
291  */
292 vtkSmartPointer<vtkPolyData> mitk::FiberBundle::GetFiberPolyData() const
293 {
294  return m_FiberPolyData;
295 }
296 
297 void mitk::FiberBundle::ColorFibersByOrientation()
298 {
299  //===== FOR WRITING A TEST ========================
300  // colorT size == tupelComponents * tupelElements
301  // compare color results
302  // to cover this code 100% also polydata needed, where colorarray already exists
303  // + one fiber with exactly 1 point
304  // + one fiber with 0 points
305  //=================================================
306 
307  vtkPoints* extrPoints = nullptr;
308  extrPoints = m_FiberPolyData->GetPoints();
309  int numOfPoints = 0;
310  if (extrPoints!=nullptr)
311  numOfPoints = extrPoints->GetNumberOfPoints();
312 
313  //colors and alpha value for each single point, RGBA = 4 components
314  unsigned char rgba[4] = {0,0,0,0};
315  int componentSize = 4;
316  m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
317  m_FiberColors->Allocate(numOfPoints * componentSize);
318  m_FiberColors->SetNumberOfComponents(componentSize);
319  m_FiberColors->SetName("FIBER_COLORS");
320 
321  int numOfFibers = m_FiberPolyData->GetNumberOfLines();
322  if (numOfFibers < 1)
323  return;
324 
325  /* extract single fibers of fiberBundle */
326  vtkCellArray* fiberList = m_FiberPolyData->GetLines();
327  fiberList->InitTraversal();
328  for (int fi=0; fi<numOfFibers; ++fi) {
329 
330  vtkIdType* idList; // contains the point id's of the line
331  vtkIdType pointsPerFiber; // number of points for current line
332  fiberList->GetNextCell(pointsPerFiber, idList);
333 
334  /* single fiber checkpoints: is number of points valid */
335  if (pointsPerFiber > 1)
336  {
337  /* operate on points of single fiber */
338  for (int i=0; i <pointsPerFiber; ++i)
339  {
340  /* process all points except starting and endpoint for calculating color value take current point, previous point and next point */
341  if (i<pointsPerFiber-1 && i > 0)
342  {
343  /* The color value of the current point is influenced by the previous point and next point. */
344  vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
345  vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]);
346  vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]);
347 
348  vnl_vector_fixed< double, 3 > diff1;
349  diff1 = currentPntvtk - nextPntvtk;
350 
351  vnl_vector_fixed< double, 3 > diff2;
352  diff2 = currentPntvtk - prevPntvtk;
353 
354  vnl_vector_fixed< double, 3 > diff;
355  diff = (diff1 - diff2) / 2.0;
356  diff.normalize();
357 
358  rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0]));
359  rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1]));
360  rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2]));
361  rgba[3] = (unsigned char) (255.0);
362  }
363  else if (i==0)
364  {
365  /* First point has no previous point, therefore only diff1 is taken */
366 
367  vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
368  vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]);
369 
370  vnl_vector_fixed< double, 3 > diff1;
371  diff1 = currentPntvtk - nextPntvtk;
372  diff1.normalize();
373 
374  rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0]));
375  rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1]));
376  rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2]));
377  rgba[3] = (unsigned char) (255.0);
378  }
379  else if (i==pointsPerFiber-1)
380  {
381  /* Last point has no next point, therefore only diff2 is taken */
382  vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]);
383  vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]);
384 
385  vnl_vector_fixed< double, 3 > diff2;
386  diff2 = currentPntvtk - prevPntvtk;
387  diff2.normalize();
388 
389  rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0]));
390  rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1]));
391  rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2]));
392  rgba[3] = (unsigned char) (255.0);
393  }
394  m_FiberColors->InsertTupleValue(idList[i], rgba);
395  }
396  }
397  else if (pointsPerFiber == 1)
398  {
399  /* a single point does not define a fiber (use vertex mechanisms instead */
400  continue;
401  }
402  else
403  {
404  MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ;
405  continue;
406  }
407  }
408  m_UpdateTime3D.Modified();
409  m_UpdateTime2D.Modified();
410 }
411 
412 void mitk::FiberBundle::ColorFibersByCurvature(bool minMaxNorm)
413 {
414  double window = 5;
415 
416  //colors and alpha value for each single point, RGBA = 4 components
417  unsigned char rgba[4] = {0,0,0,0};
418  int componentSize = 4;
419  m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
420  m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize);
421  m_FiberColors->SetNumberOfComponents(componentSize);
422  m_FiberColors->SetName("FIBER_COLORS");
423 
424  mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
425  vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
426  lookupTable->SetTableRange(0.0, 0.8);
427  lookupTable->Build();
428  mitkLookup->SetVtkLookupTable(lookupTable);
429  mitkLookup->SetType(mitk::LookupTable::JET);
430 
431  vector< double > values;
432  double min = 1;
433  double max = 0;
434  MITK_INFO << "Coloring fibers by curvature";
435  boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
436  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
437  {
438  ++disp;
439  vtkCell* cell = m_FiberPolyData->GetCell(i);
440  int numPoints = cell->GetNumberOfPoints();
441  vtkPoints* points = cell->GetPoints();
442 
443  // calculate curvatures
444  for (int j=0; j<numPoints; j++)
445  {
446  double dist = 0;
447  int c = j;
448  std::vector< vnl_vector_fixed< float, 3 > > vectors;
449  vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0);
450  while(dist<window/2 && c>1)
451  {
452  double p1[3];
453  points->GetPoint(c-1, p1);
454  double p2[3];
455  points->GetPoint(c, p2);
456 
457  vnl_vector_fixed< float, 3 > v;
458  v[0] = p2[0]-p1[0];
459  v[1] = p2[1]-p1[1];
460  v[2] = p2[2]-p1[2];
461  dist += v.magnitude();
462  v.normalize();
463  vectors.push_back(v);
464  if (c==j)
465  meanV += v;
466  c--;
467  }
468  c = j;
469  dist = 0;
470  while(dist<window/2 && c<numPoints-1)
471  {
472  double p1[3];
473  points->GetPoint(c, p1);
474  double p2[3];
475  points->GetPoint(c+1, p2);
476 
477  vnl_vector_fixed< float, 3 > v;
478  v[0] = p2[0]-p1[0];
479  v[1] = p2[1]-p1[1];
480  v[2] = p2[2]-p1[2];
481  dist += v.magnitude();
482  v.normalize();
483  vectors.push_back(v);
484  if (c==j)
485  meanV += v;
486  c++;
487  }
488  meanV.normalize();
489 
490  double dev = 0;
491  for (unsigned int c=0; c<vectors.size(); c++)
492  {
493  double angle = dot_product(meanV, vectors.at(c));
494  if (angle>1.0)
495  angle = 1.0;
496  if (angle<-1.0)
497  angle = -1.0;
498  dev += acos(angle)*180/M_PI;
499  }
500  if (vectors.size()>0)
501  dev /= vectors.size();
502 
503  dev = 1.0-dev/180.0;
504  values.push_back(dev);
505  if (dev<min)
506  min = dev;
507  if (dev>max)
508  max = dev;
509  }
510  }
511  unsigned int count = 0;
512  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
513  {
514  vtkCell* cell = m_FiberPolyData->GetCell(i);
515  int numPoints = cell->GetNumberOfPoints();
516  for (int j=0; j<numPoints; j++)
517  {
518  double color[3];
519  double dev = values.at(count);
520  if (minMaxNorm)
521  dev = (dev-min)/(max-min);
522  // double dev = values.at(count)*values.at(count);
523  lookupTable->GetColor(dev, color);
524 
525  rgba[0] = (unsigned char) (255.0 * color[0]);
526  rgba[1] = (unsigned char) (255.0 * color[1]);
527  rgba[2] = (unsigned char) (255.0 * color[2]);
528  rgba[3] = (unsigned char) (255.0);
529  m_FiberColors->InsertTupleValue(cell->GetPointId(j), rgba);
530  count++;
531  }
532  }
533  m_UpdateTime3D.Modified();
534  m_UpdateTime2D.Modified();
535 }
536 
537 void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray)
538 {
539  for(long i=0; i<m_FiberColors->GetNumberOfTuples(); i++)
540  {
541  double faValue = FAValArray->GetValue(i);
542  faValue = faValue * 255.0;
543  m_FiberColors->SetComponent(i,3, (unsigned char) faValue );
544  }
545  m_UpdateTime3D.Modified();
546  m_UpdateTime2D.Modified();
547 }
548 
549 void mitk::FiberBundle::ResetFiberOpacity()
550 {
551  for(long i=0; i<m_FiberColors->GetNumberOfTuples(); i++)
552  m_FiberColors->SetComponent(i,3, 255.0 );
553  m_UpdateTime3D.Modified();
554  m_UpdateTime2D.Modified();
555 }
556 
557 void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity)
558 {
559  mitkPixelTypeMultiplex2( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity );
560  m_UpdateTime3D.Modified();
561  m_UpdateTime2D.Modified();
562 }
563 
564 template <typename TPixel>
565 void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity)
566 {
567  m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
568  m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
569  m_FiberColors->SetNumberOfComponents(4);
570  m_FiberColors->SetName("FIBER_COLORS");
571 
572  mitk::ImagePixelReadAccessor<TPixel,3> readimage(image, image->GetVolumeData(0));
573 
574  unsigned char rgba[4] = {0,0,0,0};
575  vtkPoints* pointSet = m_FiberPolyData->GetPoints();
576 
577  mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New();
578  vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
579  lookupTable->SetTableRange(0.0, 0.8);
580  lookupTable->Build();
581  mitkLookup->SetVtkLookupTable(lookupTable);
582  mitkLookup->SetType(mitk::LookupTable::JET);
583 
584  for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
585  {
586  Point3D px;
587  px[0] = pointSet->GetPoint(i)[0];
588  px[1] = pointSet->GetPoint(i)[1];
589  px[2] = pointSet->GetPoint(i)[2];
590  double pixelValue = readimage.GetPixelByWorldCoordinates(px);
591 
592  double color[3];
593  lookupTable->GetColor(1-pixelValue, color);
594 
595  rgba[0] = (unsigned char) (255.0 * color[0]);
596  rgba[1] = (unsigned char) (255.0 * color[1]);
597  rgba[2] = (unsigned char) (255.0 * color[2]);
598  if (opacity)
599  rgba[3] = (unsigned char) (255.0 * pixelValue);
600  else
601  rgba[3] = (unsigned char) (255.0);
602  m_FiberColors->InsertTupleValue(i, rgba);
603  }
604  m_UpdateTime3D.Modified();
605  m_UpdateTime2D.Modified();
606 }
607 
608 void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha)
609 {
610  m_FiberColors = vtkSmartPointer<vtkUnsignedCharArray>::New();
611  m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4);
612  m_FiberColors->SetNumberOfComponents(4);
613  m_FiberColors->SetName("FIBER_COLORS");
614 
615  unsigned char rgba[4] = {0,0,0,0};
616  for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
617  {
618  rgba[0] = (unsigned char) r;
619  rgba[1] = (unsigned char) g;
620  rgba[2] = (unsigned char) b;
621  rgba[3] = (unsigned char) alpha;
622  m_FiberColors->InsertTupleValue(i, rgba);
623  }
624  m_UpdateTime3D.Modified();
625  m_UpdateTime2D.Modified();
626 }
627 
628 void mitk::FiberBundle::GenerateFiberIds()
629 {
630  if (m_FiberPolyData == nullptr)
631  return;
632 
633  vtkSmartPointer<vtkIdFilter> idFiberFilter = vtkSmartPointer<vtkIdFilter>::New();
634  idFiberFilter->SetInputData(m_FiberPolyData);
635  idFiberFilter->CellIdsOn();
636  // idFiberFilter->PointIdsOn(); // point id's are not needed
637  idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY);
638  idFiberFilter->FieldDataOn();
639  idFiberFilter->Update();
640 
641  m_FiberIdDataSet = idFiberFilter->GetOutput();
642 
643 }
644 
645 mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds)
646 {
647  vtkSmartPointer<vtkPolyData> polyData = m_FiberPolyData;
648  if (anyPoint)
649  {
650  float minSpacing = 1;
651  if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
652  minSpacing = mask->GetSpacing()[0];
653  else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
654  minSpacing = mask->GetSpacing()[1];
655  else
656  minSpacing = mask->GetSpacing()[2];
657 
658  mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy();
659  fibCopy->ResampleSpline(minSpacing/5);
660  polyData = fibCopy->GetFiberPolyData();
661  }
662  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
663  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
664 
665  MITK_INFO << "Extracting fibers";
666  boost::progress_display disp(m_NumFibers);
667  for (int i=0; i<m_NumFibers; i++)
668  {
669  ++disp;
670 
671  vtkCell* cell = polyData->GetCell(i);
672  int numPoints = cell->GetNumberOfPoints();
673  vtkPoints* points = cell->GetPoints();
674 
675  vtkCell* cellOriginal = m_FiberPolyData->GetCell(i);
676  int numPointsOriginal = cellOriginal->GetNumberOfPoints();
677  vtkPoints* pointsOriginal = cellOriginal->GetPoints();
678 
679  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
680 
681  if (numPoints>1 && numPointsOriginal)
682  {
683  if (anyPoint)
684  {
685  if (!invert)
686  {
687  for (int j=0; j<numPoints; j++)
688  {
689  double* p = points->GetPoint(j);
690 
691  itk::Point<float, 3> itkP;
692  itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
693  itk::Index<3> idx;
694  mask->TransformPhysicalPointToIndex(itkP, idx);
695 
696  if ( mask->GetPixel(idx)>0 && mask->GetLargestPossibleRegion().IsInside(idx) )
697  {
698  for (int k=0; k<numPointsOriginal; k++)
699  {
700  double* p = pointsOriginal->GetPoint(k);
701  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
702  container->GetPointIds()->InsertNextId(id);
703  }
704  break;
705  }
706  }
707  }
708  else
709  {
710  bool includeFiber = true;
711  for (int j=0; j<numPoints; j++)
712  {
713  double* p = points->GetPoint(j);
714 
715  itk::Point<float, 3> itkP;
716  itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
717  itk::Index<3> idx;
718  mask->TransformPhysicalPointToIndex(itkP, idx);
719 
720  if ( mask->GetPixel(idx)>0 && mask->GetLargestPossibleRegion().IsInside(idx) )
721  {
722  includeFiber = false;
723  break;
724  }
725  }
726  if (includeFiber)
727  {
728 
729  for (int k=0; k<numPointsOriginal; k++)
730  {
731  double* p = pointsOriginal->GetPoint(k);
732  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
733  container->GetPointIds()->InsertNextId(id);
734  }
735  }
736  }
737  }
738  else
739  {
740  double* start = pointsOriginal->GetPoint(0);
741  itk::Point<float, 3> itkStart;
742  itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2];
743  itk::Index<3> idxStart;
744  mask->TransformPhysicalPointToIndex(itkStart, idxStart);
745 
746  double* end = pointsOriginal->GetPoint(numPointsOriginal-1);
747  itk::Point<float, 3> itkEnd;
748  itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2];
749  itk::Index<3> idxEnd;
750  mask->TransformPhysicalPointToIndex(itkEnd, idxEnd);
751 
752  if (invert)
753  {
754  if (bothEnds)
755  {
756  if ( !mask->GetPixel(idxStart)>0 && !mask->GetPixel(idxEnd)>0 )
757  {
758  for (int j=0; j<numPointsOriginal; j++)
759  {
760  double* p = pointsOriginal->GetPoint(j);
761  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
762  container->GetPointIds()->InsertNextId(id);
763  }
764  }
765  }
766  else if ( !mask->GetPixel(idxStart)>0 || !mask->GetPixel(idxEnd)>0 )
767  {
768  for (int j=0; j<numPointsOriginal; j++)
769  {
770  double* p = pointsOriginal->GetPoint(j);
771  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
772  container->GetPointIds()->InsertNextId(id);
773  }
774  }
775  }
776  else
777  {
778  if (bothEnds)
779  {
780  if ( mask->GetPixel(idxStart)>0 && mask->GetPixel(idxEnd)>0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) )
781  {
782  for (int j=0; j<numPointsOriginal; j++)
783  {
784  double* p = pointsOriginal->GetPoint(j);
785  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
786  container->GetPointIds()->InsertNextId(id);
787  }
788  }
789  }
790  else if ( (mask->GetPixel(idxStart)>0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd)>0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) )
791  {
792  for (int j=0; j<numPointsOriginal; j++)
793  {
794  double* p = pointsOriginal->GetPoint(j);
795  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
796  container->GetPointIds()->InsertNextId(id);
797  }
798  }
799  }
800  }
801  }
802 
803  vtkNewCells->InsertNextCell(container);
804  }
805 
806  if (vtkNewCells->GetNumberOfCells()<=0)
807  return nullptr;
808 
809  vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
810  newPolyData->SetPoints(vtkNewPoints);
811  newPolyData->SetLines(vtkNewCells);
812  return mitk::FiberBundle::New(newPolyData);
813 }
814 
815 mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert)
816 {
817  float minSpacing = 1;
818  if(mask->GetSpacing()[0]<mask->GetSpacing()[1] && mask->GetSpacing()[0]<mask->GetSpacing()[2])
819  minSpacing = mask->GetSpacing()[0];
820  else if (mask->GetSpacing()[1] < mask->GetSpacing()[2])
821  minSpacing = mask->GetSpacing()[1];
822  else
823  minSpacing = mask->GetSpacing()[2];
824 
825  mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy();
826  fibCopy->ResampleSpline(minSpacing/10);
827  vtkSmartPointer<vtkPolyData> polyData =fibCopy->GetFiberPolyData();
828 
829  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
830  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
831 
832  MITK_INFO << "Cutting fibers";
833  boost::progress_display disp(m_NumFibers);
834  for (int i=0; i<m_NumFibers; i++)
835  {
836  ++disp;
837 
838  vtkCell* cell = polyData->GetCell(i);
839  int numPoints = cell->GetNumberOfPoints();
840  vtkPoints* points = cell->GetPoints();
841 
842  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
843  if (numPoints>1)
844  {
845  int newNumPoints = 0;
846  for (int j=0; j<numPoints; j++)
847  {
848  double* p = points->GetPoint(j);
849 
850  itk::Point<float, 3> itkP;
851  itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2];
852  itk::Index<3> idx;
853  mask->TransformPhysicalPointToIndex(itkP, idx);
854 
855  if ( mask->GetPixel(idx)>0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert )
856  {
857  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
858  container->GetPointIds()->InsertNextId(id);
859  newNumPoints++;
860  }
861  else if ( (mask->GetPixel(idx)<=0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert )
862  {
863  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
864  container->GetPointIds()->InsertNextId(id);
865  newNumPoints++;
866  }
867  else if (newNumPoints>0)
868  {
869  vtkNewCells->InsertNextCell(container);
870 
871  newNumPoints = 0;
872  container = vtkSmartPointer<vtkPolyLine>::New();
873  }
874  }
875 
876  if (newNumPoints>0)
877  vtkNewCells->InsertNextCell(container);
878  }
879 
880  }
881 
882  if (vtkNewCells->GetNumberOfCells()<=0)
883  return nullptr;
884 
885  vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
886  newPolyData->SetPoints(vtkNewPoints);
887  newPolyData->SetLines(vtkNewCells);
888  mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData);
889  newFib->Compress(0.1);
890  return newFib;
891 }
892 
893 mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage)
894 {
895  if (roi==nullptr || !(dynamic_cast<PlanarFigure*>(roi->GetData()) || dynamic_cast<PlanarFigureComposite*>(roi->GetData())) )
896  return nullptr;
897 
898  std::vector<long> tmp = ExtractFiberIdSubset(roi, storage);
899 
900  if (tmp.size()<=0)
901  return mitk::FiberBundle::New();
902  vtkSmartPointer<vtkPolyData> pTmp = GeneratePolyDataByIds(tmp);
903  return mitk::FiberBundle::New(pTmp);
904 }
905 
906 std::vector<long> mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage)
907 {
908  std::vector<long> result;
909  if (roi==nullptr || roi->GetData()==nullptr)
910  return result;
911 
912  mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast<mitk::PlanarFigureComposite*>(roi->GetData());
913  if (!pfc.IsNull()) // handle composite
914  {
915  DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi);
916  if (children->size()==0)
917  return result;
918 
919  switch (pfc->getOperationType())
920  {
921  case 0: // AND
922  {
923  MITK_INFO << "AND";
924  result = this->ExtractFiberIdSubset(children->ElementAt(0), storage);
925  std::vector<long>::iterator it;
926  for (unsigned int i=1; i<children->Size(); ++i)
927  {
928  std::vector<long> inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage);
929 
930  std::vector<long> rest(std::min(result.size(),inRoi.size()));
931  it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
932  rest.resize( it - rest.begin() );
933  result = rest;
934  }
935  break;
936  }
937  case 1: // OR
938  {
939  MITK_INFO << "OR";
940  result = ExtractFiberIdSubset(children->ElementAt(0), storage);
941  std::vector<long>::iterator it;
942  for (unsigned int i=1; i<children->Size(); ++i)
943  {
944  it = result.end();
945  std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
946  result.insert(it, inRoi.begin(), inRoi.end());
947  }
948 
949  // remove duplicates
950  sort(result.begin(), result.end());
951  it = unique(result.begin(), result.end());
952  result.resize( it - result.begin() );
953  break;
954  }
955  case 2: // NOT
956  {
957  MITK_INFO << "NOT";
958  for(long i=0; i<this->GetNumFibers(); i++)
959  result.push_back(i);
960 
961  std::vector<long>::iterator it;
962  for (long i=0; i<children->Size(); ++i)
963  {
964  std::vector<long> inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage);
965 
966  std::vector<long> rest(result.size()-inRoi.size());
967  it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() );
968  rest.resize( it - rest.begin() );
969  result = rest;
970  }
971  break;
972  }
973  }
974  }
975  else if ( dynamic_cast<mitk::PlanarFigure*>(roi->GetData()) ) // actual extraction
976  {
977  if ( dynamic_cast<mitk::PlanarPolygon*>(roi->GetData()) )
978  {
979  mitk::PlanarFigure::Pointer planarPoly = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
980 
981  //create vtkPolygon using controlpoints from planarFigure polygon
982  vtkSmartPointer<vtkPolygon> polygonVtk = vtkSmartPointer<vtkPolygon>::New();
983  for (unsigned int i=0; i<planarPoly->GetNumberOfControlPoints(); ++i)
984  {
985  itk::Point<double,3> p = planarPoly->GetWorldControlPoint(i);
986  vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] );
987  polygonVtk->GetPointIds()->InsertNextId(id);
988  }
989 
990  MITK_INFO << "Extracting with polygon";
991  boost::progress_display disp(m_NumFibers);
992  for (int i=0; i<m_NumFibers; i++)
993  {
994  ++disp ;
995  vtkCell* cell = m_FiberPolyData->GetCell(i);
996  int numPoints = cell->GetNumberOfPoints();
997  vtkPoints* points = cell->GetPoints();
998 
999  for (int j=0; j<numPoints-1; j++)
1000  {
1001  // Inputs
1002  double p1[3] = {0,0,0};
1003  points->GetPoint(j, p1);
1004  double p2[3] = {0,0,0};
1005  points->GetPoint(j+1, p2);
1006  double tolerance = 0.001;
1007 
1008  // Outputs
1009  double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
1010  double x[3] = {0,0,0}; // The coordinate of the intersection
1011  double pcoords[3] = {0,0,0};
1012  int subId = 0;
1013 
1014  int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId);
1015  if (iD!=0)
1016  {
1017  result.push_back(i);
1018  break;
1019  }
1020  }
1021  }
1022  }
1023  else if ( dynamic_cast<mitk::PlanarCircle*>(roi->GetData()) )
1024  {
1025  mitk::PlanarFigure::Pointer planarFigure = dynamic_cast<mitk::PlanarFigure*>(roi->GetData());
1026  Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal();
1027  planeNormal.Normalize();
1028 
1029  //calculate circle radius
1030  mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint
1031  mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint
1032 
1033  double radius = V1w.EuclideanDistanceTo(V2w);
1034  radius *= radius;
1035 
1036  MITK_INFO << "Extracting with circle";
1037  boost::progress_display disp(m_NumFibers);
1038  for (int i=0; i<m_NumFibers; i++)
1039  {
1040  ++disp ;
1041  vtkCell* cell = m_FiberPolyData->GetCell(i);
1042  int numPoints = cell->GetNumberOfPoints();
1043  vtkPoints* points = cell->GetPoints();
1044 
1045  for (int j=0; j<numPoints-1; j++)
1046  {
1047  // Inputs
1048  double p1[3] = {0,0,0};
1049  points->GetPoint(j, p1);
1050  double p2[3] = {0,0,0};
1051  points->GetPoint(j+1, p2);
1052 
1053  // Outputs
1054  double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2))
1055  double x[3] = {0,0,0}; // The coordinate of the intersection
1056 
1057  int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x);
1058 
1059  if (iD!=0)
1060  {
1061  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]);
1062  if( dist <= radius)
1063  {
1064  result.push_back(i);
1065  break;
1066  }
1067  }
1068  }
1069  }
1070  }
1071  return result;
1072  }
1073 
1074  return result;
1075 }
1076 
1077 void mitk::FiberBundle::UpdateFiberGeometry()
1078 {
1079  vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
1080  cleaner->SetInputData(m_FiberPolyData);
1081  cleaner->PointMergingOff();
1082  cleaner->Update();
1083  m_FiberPolyData = cleaner->GetOutput();
1084 
1085  m_FiberLengths.clear();
1086  m_MeanFiberLength = 0;
1087  m_MedianFiberLength = 0;
1088  m_LengthStDev = 0;
1089  m_NumFibers = m_FiberPolyData->GetNumberOfCells();
1090 
1091  if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints())
1092  this->ColorFibersByOrientation();
1093 
1094  if (m_FiberWeights->GetSize()!=m_NumFibers)
1095  {
1096  m_FiberWeights = vtkSmartPointer<vtkFloatArray>::New();
1097  m_FiberWeights->SetName("FIBER_WEIGHTS");
1098  m_FiberWeights->SetNumberOfValues(m_NumFibers);
1099  this->SetFiberWeights(1);
1100  }
1101 
1102  if (m_NumFibers<=0) // no fibers present; apply default geometry
1103  {
1104  m_MinFiberLength = 0;
1105  m_MaxFiberLength = 0;
1106  mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
1107  geometry->SetImageGeometry(false);
1108  float b[] = {0, 1, 0, 1, 0, 1};
1109  geometry->SetFloatBounds(b);
1110  SetGeometry(geometry);
1111  return;
1112  }
1113  double b[6];
1114  m_FiberPolyData->GetBounds(b);
1115 
1116  // calculate statistics
1117  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
1118  {
1119  vtkCell* cell = m_FiberPolyData->GetCell(i);
1120  int p = cell->GetNumberOfPoints();
1121  vtkPoints* points = cell->GetPoints();
1122  float length = 0;
1123  for (int j=0; j<p-1; j++)
1124  {
1125  double p1[3];
1126  points->GetPoint(j, p1);
1127  double p2[3];
1128  points->GetPoint(j+1, p2);
1129 
1130  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]));
1131  length += dist;
1132  }
1133  m_FiberLengths.push_back(length);
1134  m_MeanFiberLength += length;
1135  if (i==0)
1136  {
1137  m_MinFiberLength = length;
1138  m_MaxFiberLength = length;
1139  }
1140  else
1141  {
1142  if (length<m_MinFiberLength)
1143  m_MinFiberLength = length;
1144  if (length>m_MaxFiberLength)
1145  m_MaxFiberLength = length;
1146  }
1147  }
1148  m_MeanFiberLength /= m_NumFibers;
1149 
1150  std::vector< float > sortedLengths = m_FiberLengths;
1151  std::sort(sortedLengths.begin(), sortedLengths.end());
1152  for (int i=0; i<m_NumFibers; i++)
1153  m_LengthStDev += (m_MeanFiberLength-sortedLengths.at(i))*(m_MeanFiberLength-sortedLengths.at(i));
1154  if (m_NumFibers>1)
1155  m_LengthStDev /= (m_NumFibers-1);
1156  else
1157  m_LengthStDev = 0;
1158  m_LengthStDev = std::sqrt(m_LengthStDev);
1159  m_MedianFiberLength = sortedLengths.at(m_NumFibers/2);
1160 
1161  mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
1162  geometry->SetFloatBounds(b);
1163  this->SetGeometry(geometry);
1164 
1165  m_UpdateTime3D.Modified();
1166  m_UpdateTime2D.Modified();
1167 }
1168 
1169 float mitk::FiberBundle::GetFiberWeight(unsigned int fiber)
1170 {
1171  return m_FiberWeights->GetValue(fiber);
1172 }
1173 
1174 void mitk::FiberBundle::SetFiberWeights(float newWeight)
1175 {
1176  for (int i=0; i<m_FiberWeights->GetSize(); i++)
1177  m_FiberWeights->SetValue(i, newWeight);
1178 }
1179 
1180 void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer<vtkFloatArray> weights)
1181 {
1182  if (m_NumFibers!=weights->GetSize())
1183  {
1184  MITK_INFO << "Weights array not equal to number of fibers!";
1185  return;
1186  }
1187 
1188  for (int i=0; i<weights->GetSize(); i++)
1189  m_FiberWeights->SetValue(i, weights->GetValue(i));
1190 
1191  m_FiberWeights->SetName("FIBER_WEIGHTS");
1192 }
1193 
1194 void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight)
1195 {
1196  m_FiberWeights->SetValue(fiber, weight);
1197 }
1198 
1199 void mitk::FiberBundle::SetFiberColors(vtkSmartPointer<vtkUnsignedCharArray> fiberColors)
1200 {
1201  for(long i=0; i<m_FiberPolyData->GetNumberOfPoints(); ++i)
1202  {
1203  unsigned char source[4] = {0,0,0,0};
1204  fiberColors->GetTupleValue(i, source);
1205 
1206  unsigned char target[4] = {0,0,0,0};
1207  target[0] = source[0];
1208  target[1] = source[1];
1209  target[2] = source[2];
1210  target[3] = source[3];
1211  m_FiberColors->InsertTupleValue(i, target);
1212  }
1213  m_UpdateTime3D.Modified();
1214  m_UpdateTime2D.Modified();
1215 }
1216 
1217 itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz)
1218 {
1219  rx = rx*M_PI/180;
1220  ry = ry*M_PI/180;
1221  rz = rz*M_PI/180;
1222 
1223  itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity();
1224  rotX[1][1] = cos(rx);
1225  rotX[2][2] = rotX[1][1];
1226  rotX[1][2] = -sin(rx);
1227  rotX[2][1] = -rotX[1][2];
1228 
1229  itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity();
1230  rotY[0][0] = cos(ry);
1231  rotY[2][2] = rotY[0][0];
1232  rotY[0][2] = sin(ry);
1233  rotY[2][0] = -rotY[0][2];
1234 
1235  itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity();
1236  rotZ[0][0] = cos(rz);
1237  rotZ[1][1] = rotZ[0][0];
1238  rotZ[0][1] = -sin(rz);
1239  rotZ[1][0] = -rotZ[0][1];
1240 
1241  itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX;
1242 
1243  m = rot*m;
1244 
1245  return m;
1246 }
1247 
1248 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)
1249 {
1250  rx = rx*M_PI/180;
1251  ry = ry*M_PI/180;
1252  rz = rz*M_PI/180;
1253 
1254  vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
1255  rotX[1][1] = cos(rx);
1256  rotX[2][2] = rotX[1][1];
1257  rotX[1][2] = -sin(rx);
1258  rotX[2][1] = -rotX[1][2];
1259 
1260  vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
1261  rotY[0][0] = cos(ry);
1262  rotY[2][2] = rotY[0][0];
1263  rotY[0][2] = sin(ry);
1264  rotY[2][0] = -rotY[0][2];
1265 
1266  vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
1267  rotZ[0][0] = cos(rz);
1268  rotZ[1][1] = rotZ[0][0];
1269  rotZ[0][1] = -sin(rz);
1270  rotZ[1][0] = -rotZ[0][1];
1271 
1272  vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX;
1273 
1274  mitk::BaseGeometry::Pointer geom = this->GetGeometry();
1275  mitk::Point3D center = geom->GetCenter();
1276 
1277  point[0] -= center[0];
1278  point[1] -= center[1];
1279  point[2] -= center[2];
1280  point = rot*point;
1281  point[0] += center[0]+tx;
1282  point[1] += center[1]+ty;
1283  point[2] += center[2]+tz;
1284  itk::Point<float, 3> out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2];
1285  return out;
1286 }
1287 
1288 void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz)
1289 {
1290  rx = rx*M_PI/180;
1291  ry = ry*M_PI/180;
1292  rz = rz*M_PI/180;
1293 
1294  vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
1295  rotX[1][1] = cos(rx);
1296  rotX[2][2] = rotX[1][1];
1297  rotX[1][2] = -sin(rx);
1298  rotX[2][1] = -rotX[1][2];
1299 
1300  vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
1301  rotY[0][0] = cos(ry);
1302  rotY[2][2] = rotY[0][0];
1303  rotY[0][2] = sin(ry);
1304  rotY[2][0] = -rotY[0][2];
1305 
1306  vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
1307  rotZ[0][0] = cos(rz);
1308  rotZ[1][1] = rotZ[0][0];
1309  rotZ[0][1] = -sin(rz);
1310  rotZ[1][0] = -rotZ[0][1];
1311 
1312  vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX;
1313 
1314  mitk::BaseGeometry::Pointer geom = this->GetGeometry();
1315  mitk::Point3D center = geom->GetCenter();
1316 
1317  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1318  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1319 
1320  for (int i=0; i<m_NumFibers; i++)
1321  {
1322  vtkCell* cell = m_FiberPolyData->GetCell(i);
1323  int numPoints = cell->GetNumberOfPoints();
1324  vtkPoints* points = cell->GetPoints();
1325 
1326  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1327  for (int j=0; j<numPoints; j++)
1328  {
1329  double* p = points->GetPoint(j);
1330  vnl_vector_fixed< double, 3 > dir;
1331  dir[0] = p[0]-center[0];
1332  dir[1] = p[1]-center[1];
1333  dir[2] = p[2]-center[2];
1334  dir = rot*dir;
1335  dir[0] += center[0]+tx;
1336  dir[1] += center[1]+ty;
1337  dir[2] += center[2]+tz;
1338  vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block());
1339  container->GetPointIds()->InsertNextId(id);
1340  }
1341  vtkNewCells->InsertNextCell(container);
1342  }
1343 
1344  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1345  m_FiberPolyData->SetPoints(vtkNewPoints);
1346  m_FiberPolyData->SetLines(vtkNewCells);
1347  this->SetFiberPolyData(m_FiberPolyData, true);
1348 }
1349 
1350 void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z)
1351 {
1352  x = x*M_PI/180;
1353  y = y*M_PI/180;
1354  z = z*M_PI/180;
1355 
1356  vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity();
1357  rotX[1][1] = cos(x);
1358  rotX[2][2] = rotX[1][1];
1359  rotX[1][2] = -sin(x);
1360  rotX[2][1] = -rotX[1][2];
1361 
1362  vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity();
1363  rotY[0][0] = cos(y);
1364  rotY[2][2] = rotY[0][0];
1365  rotY[0][2] = sin(y);
1366  rotY[2][0] = -rotY[0][2];
1367 
1368  vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
1369  rotZ[0][0] = cos(z);
1370  rotZ[1][1] = rotZ[0][0];
1371  rotZ[0][1] = -sin(z);
1372  rotZ[1][0] = -rotZ[0][1];
1373 
1374  mitk::BaseGeometry::Pointer geom = this->GetGeometry();
1375  mitk::Point3D center = geom->GetCenter();
1376 
1377  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1378  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1379 
1380  for (int i=0; i<m_NumFibers; i++)
1381  {
1382  vtkCell* cell = m_FiberPolyData->GetCell(i);
1383  int numPoints = cell->GetNumberOfPoints();
1384  vtkPoints* points = cell->GetPoints();
1385 
1386  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1387  for (int j=0; j<numPoints; j++)
1388  {
1389  double* p = points->GetPoint(j);
1390  vnl_vector_fixed< double, 3 > dir;
1391  dir[0] = p[0]-center[0];
1392  dir[1] = p[1]-center[1];
1393  dir[2] = p[2]-center[2];
1394  dir = rotZ*rotY*rotX*dir;
1395  dir[0] += center[0];
1396  dir[1] += center[1];
1397  dir[2] += center[2];
1398  vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block());
1399  container->GetPointIds()->InsertNextId(id);
1400  }
1401  vtkNewCells->InsertNextCell(container);
1402  }
1403 
1404  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1405  m_FiberPolyData->SetPoints(vtkNewPoints);
1406  m_FiberPolyData->SetLines(vtkNewCells);
1407  this->SetFiberPolyData(m_FiberPolyData, true);
1408 }
1409 
1410 void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter)
1411 {
1412  MITK_INFO << "Scaling fibers";
1413  boost::progress_display disp(m_NumFibers);
1414 
1415  mitk::BaseGeometry* geom = this->GetGeometry();
1416  mitk::Point3D c = geom->GetCenter();
1417 
1418  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1419  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1420 
1421  for (int i=0; i<m_NumFibers; i++)
1422  {
1423  ++disp ;
1424  vtkCell* cell = m_FiberPolyData->GetCell(i);
1425  int numPoints = cell->GetNumberOfPoints();
1426  vtkPoints* points = cell->GetPoints();
1427 
1428  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1429  for (int j=0; j<numPoints; j++)
1430  {
1431  double* p = points->GetPoint(j);
1432  if (subtractCenter)
1433  {
1434  p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2];
1435  }
1436  p[0] *= x;
1437  p[1] *= y;
1438  p[2] *= z;
1439  if (subtractCenter)
1440  {
1441  p[0] += c[0]; p[1] += c[1]; p[2] += c[2];
1442  }
1443  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
1444  container->GetPointIds()->InsertNextId(id);
1445  }
1446  vtkNewCells->InsertNextCell(container);
1447  }
1448 
1449  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1450  m_FiberPolyData->SetPoints(vtkNewPoints);
1451  m_FiberPolyData->SetLines(vtkNewCells);
1452  this->SetFiberPolyData(m_FiberPolyData, true);
1453 }
1454 
1455 void mitk::FiberBundle::TranslateFibers(double x, double y, double z)
1456 {
1457  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1458  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1459 
1460  for (int i=0; i<m_NumFibers; i++)
1461  {
1462  vtkCell* cell = m_FiberPolyData->GetCell(i);
1463  int numPoints = cell->GetNumberOfPoints();
1464  vtkPoints* points = cell->GetPoints();
1465 
1466  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1467  for (int j=0; j<numPoints; j++)
1468  {
1469  double* p = points->GetPoint(j);
1470  p[0] += x;
1471  p[1] += y;
1472  p[2] += z;
1473  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
1474  container->GetPointIds()->InsertNextId(id);
1475  }
1476  vtkNewCells->InsertNextCell(container);
1477  }
1478 
1479  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1480  m_FiberPolyData->SetPoints(vtkNewPoints);
1481  m_FiberPolyData->SetLines(vtkNewCells);
1482  this->SetFiberPolyData(m_FiberPolyData, true);
1483 }
1484 
1485 void mitk::FiberBundle::MirrorFibers(unsigned int axis)
1486 {
1487  if (axis>2)
1488  return;
1489 
1490  MITK_INFO << "Mirroring fibers";
1491  boost::progress_display disp(m_NumFibers);
1492 
1493  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1494  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1495 
1496  for (int i=0; i<m_NumFibers; i++)
1497  {
1498  ++disp;
1499  vtkCell* cell = m_FiberPolyData->GetCell(i);
1500  int numPoints = cell->GetNumberOfPoints();
1501  vtkPoints* points = cell->GetPoints();
1502 
1503  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1504  for (int j=0; j<numPoints; j++)
1505  {
1506  double* p = points->GetPoint(j);
1507  p[axis] = -p[axis];
1508  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
1509  container->GetPointIds()->InsertNextId(id);
1510  }
1511  vtkNewCells->InsertNextCell(container);
1512  }
1513 
1514  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1515  m_FiberPolyData->SetPoints(vtkNewPoints);
1516  m_FiberPolyData->SetLines(vtkNewCells);
1517  this->SetFiberPolyData(m_FiberPolyData, true);
1518 }
1519 
1520 void mitk::FiberBundle::RemoveDir(vnl_vector_fixed<double,3> dir, double threshold)
1521 {
1522  dir.normalize();
1523  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1524  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1525 
1526  boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
1527  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
1528  {
1529  ++disp ;
1530  vtkCell* cell = m_FiberPolyData->GetCell(i);
1531  int numPoints = cell->GetNumberOfPoints();
1532  vtkPoints* points = cell->GetPoints();
1533 
1534  // calculate curvatures
1535  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1536  bool discard = false;
1537  for (int j=0; j<numPoints-1; j++)
1538  {
1539  double p1[3];
1540  points->GetPoint(j, p1);
1541  double p2[3];
1542  points->GetPoint(j+1, p2);
1543 
1544  vnl_vector_fixed< double, 3 > v1;
1545  v1[0] = p2[0]-p1[0];
1546  v1[1] = p2[1]-p1[1];
1547  v1[2] = p2[2]-p1[2];
1548  if (v1.magnitude()>0.001)
1549  {
1550  v1.normalize();
1551 
1552  if (fabs(dot_product(v1,dir))>threshold)
1553  {
1554  discard = true;
1555  break;
1556  }
1557  }
1558  }
1559  if (!discard)
1560  {
1561  for (int j=0; j<numPoints; j++)
1562  {
1563  double p1[3];
1564  points->GetPoint(j, p1);
1565 
1566  vtkIdType id = vtkNewPoints->InsertNextPoint(p1);
1567  container->GetPointIds()->InsertNextId(id);
1568  }
1569  vtkNewCells->InsertNextCell(container);
1570  }
1571  }
1572 
1573  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1574  m_FiberPolyData->SetPoints(vtkNewPoints);
1575  m_FiberPolyData->SetLines(vtkNewCells);
1576 
1577  this->SetFiberPolyData(m_FiberPolyData, true);
1578 
1579  // UpdateColorCoding();
1580  // UpdateFiberGeometry();
1581 }
1582 
1583 bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers)
1584 {
1585  if (minRadius<0)
1586  return true;
1587 
1588  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1589  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1590 
1591  MITK_INFO << "Applying curvature threshold";
1592  boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
1593  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
1594  {
1595  ++disp ;
1596  vtkCell* cell = m_FiberPolyData->GetCell(i);
1597  int numPoints = cell->GetNumberOfPoints();
1598  vtkPoints* points = cell->GetPoints();
1599 
1600  // calculate curvatures
1601  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1602  for (int j=0; j<numPoints-2; j++)
1603  {
1604  double p1[3];
1605  points->GetPoint(j, p1);
1606  double p2[3];
1607  points->GetPoint(j+1, p2);
1608  double p3[3];
1609  points->GetPoint(j+2, p3);
1610 
1611  vnl_vector_fixed< float, 3 > v1, v2, v3;
1612 
1613  v1[0] = p2[0]-p1[0];
1614  v1[1] = p2[1]-p1[1];
1615  v1[2] = p2[2]-p1[2];
1616 
1617  v2[0] = p3[0]-p2[0];
1618  v2[1] = p3[1]-p2[1];
1619  v2[2] = p3[2]-p2[2];
1620 
1621  v3[0] = p1[0]-p3[0];
1622  v3[1] = p1[1]-p3[1];
1623  v3[2] = p1[2]-p3[2];
1624 
1625  float a = v1.magnitude();
1626  float b = v2.magnitude();
1627  float c = v3.magnitude();
1628  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)
1629 
1630  vtkIdType id = vtkNewPoints->InsertNextPoint(p1);
1631  container->GetPointIds()->InsertNextId(id);
1632 
1633  if (deleteFibers && r<minRadius)
1634  break;
1635 
1636  if (r<minRadius)
1637  {
1638  j += 2;
1639  vtkNewCells->InsertNextCell(container);
1640  container = vtkSmartPointer<vtkPolyLine>::New();
1641  }
1642  else if (j==numPoints-3)
1643  {
1644  id = vtkNewPoints->InsertNextPoint(p2);
1645  container->GetPointIds()->InsertNextId(id);
1646  id = vtkNewPoints->InsertNextPoint(p3);
1647  container->GetPointIds()->InsertNextId(id);
1648  vtkNewCells->InsertNextCell(container);
1649  }
1650  }
1651  }
1652 
1653  if (vtkNewCells->GetNumberOfCells()<=0)
1654  return false;
1655 
1656  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1657  m_FiberPolyData->SetPoints(vtkNewPoints);
1658  m_FiberPolyData->SetLines(vtkNewCells);
1659  this->SetFiberPolyData(m_FiberPolyData, true);
1660  return true;
1661 }
1662 
1663 bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM)
1664 {
1665  MITK_INFO << "Removing short fibers";
1666  if (lengthInMM<=0 || lengthInMM<m_MinFiberLength)
1667  {
1668  MITK_INFO << "No fibers shorter than " << lengthInMM << " mm found!";
1669  return true;
1670  }
1671 
1672  if (lengthInMM>m_MaxFiberLength) // can't remove all fibers
1673  {
1674  MITK_WARN << "Process aborted. No fibers would be left!";
1675  return false;
1676  }
1677 
1678  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1679  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1680  float min = m_MaxFiberLength;
1681 
1682  boost::progress_display disp(m_NumFibers);
1683  for (int i=0; i<m_NumFibers; i++)
1684  {
1685  ++disp;
1686  vtkCell* cell = m_FiberPolyData->GetCell(i);
1687  int numPoints = cell->GetNumberOfPoints();
1688  vtkPoints* points = cell->GetPoints();
1689 
1690  if (m_FiberLengths.at(i)>=lengthInMM)
1691  {
1692  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1693  for (int j=0; j<numPoints; j++)
1694  {
1695  double* p = points->GetPoint(j);
1696  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
1697  container->GetPointIds()->InsertNextId(id);
1698  }
1699  vtkNewCells->InsertNextCell(container);
1700  if (m_FiberLengths.at(i)<min)
1701  min = m_FiberLengths.at(i);
1702  }
1703  }
1704 
1705  if (vtkNewCells->GetNumberOfCells()<=0)
1706  return false;
1707 
1708  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1709  m_FiberPolyData->SetPoints(vtkNewPoints);
1710  m_FiberPolyData->SetLines(vtkNewCells);
1711  this->SetFiberPolyData(m_FiberPolyData, true);
1712  return true;
1713 }
1714 
1715 bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM)
1716 {
1717  if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength)
1718  return true;
1719 
1720  if (lengthInMM<m_MinFiberLength) // can't remove all fibers
1721  return false;
1722 
1723  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1724  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1725 
1726  MITK_INFO << "Removing long fibers";
1727  boost::progress_display disp(m_NumFibers);
1728  for (int i=0; i<m_NumFibers; i++)
1729  {
1730  ++disp;
1731  vtkCell* cell = m_FiberPolyData->GetCell(i);
1732  int numPoints = cell->GetNumberOfPoints();
1733  vtkPoints* points = cell->GetPoints();
1734 
1735  if (m_FiberLengths.at(i)<=lengthInMM)
1736  {
1737  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1738  for (int j=0; j<numPoints; j++)
1739  {
1740  double* p = points->GetPoint(j);
1741  vtkIdType id = vtkNewPoints->InsertNextPoint(p);
1742  container->GetPointIds()->InsertNextId(id);
1743  }
1744  vtkNewCells->InsertNextCell(container);
1745  }
1746  }
1747 
1748  if (vtkNewCells->GetNumberOfCells()<=0)
1749  return false;
1750 
1751  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1752  m_FiberPolyData->SetPoints(vtkNewPoints);
1753  m_FiberPolyData->SetLines(vtkNewCells);
1754  this->SetFiberPolyData(m_FiberPolyData, true);
1755  return true;
1756 }
1757 
1758 void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias )
1759 {
1760  if (pointDistance<=0)
1761  return;
1762 
1763  vtkSmartPointer<vtkPoints> vtkSmoothPoints = vtkSmartPointer<vtkPoints>::New(); //in smoothpoints the interpolated points representing a fiber are stored.
1764 
1765  //in vtkcells all polylines are stored, actually all id's of them are stored
1766  vtkSmartPointer<vtkCellArray> vtkSmoothCells = vtkSmartPointer<vtkCellArray>::New(); //cellcontainer for smoothed lines
1767  vtkIdType pointHelperCnt = 0;
1768 
1769  MITK_INFO << "Smoothing fibers";
1770  boost::progress_display disp(m_NumFibers);
1771  for (int i=0; i<m_NumFibers; i++)
1772  {
1773  ++disp;
1774  vtkCell* cell = m_FiberPolyData->GetCell(i);
1775  int numPoints = cell->GetNumberOfPoints();
1776  vtkPoints* points = cell->GetPoints();
1777 
1778  vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
1779  for (int j=0; j<numPoints; j++)
1780  newPoints->InsertNextPoint(points->GetPoint(j));
1781 
1782  float length = m_FiberLengths.at(i);
1783  int sampling = std::ceil(length/pointDistance);
1784 
1785  vtkSmartPointer<vtkKochanekSpline> xSpline = vtkSmartPointer<vtkKochanekSpline>::New();
1786  vtkSmartPointer<vtkKochanekSpline> ySpline = vtkSmartPointer<vtkKochanekSpline>::New();
1787  vtkSmartPointer<vtkKochanekSpline> zSpline = vtkSmartPointer<vtkKochanekSpline>::New();
1788  xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity);
1789  ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity);
1790  zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity);
1791 
1792  vtkSmartPointer<vtkParametricSpline> spline = vtkSmartPointer<vtkParametricSpline>::New();
1793  spline->SetXSpline(xSpline);
1794  spline->SetYSpline(ySpline);
1795  spline->SetZSpline(zSpline);
1796  spline->SetPoints(newPoints);
1797 
1798  vtkSmartPointer<vtkParametricFunctionSource> functionSource = vtkSmartPointer<vtkParametricFunctionSource>::New();
1799  functionSource->SetParametricFunction(spline);
1800  functionSource->SetUResolution(sampling);
1801  functionSource->SetVResolution(sampling);
1802  functionSource->SetWResolution(sampling);
1803  functionSource->Update();
1804 
1805  vtkPolyData* outputFunction = functionSource->GetOutput();
1806  vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber
1807 
1808  vtkSmartPointer<vtkPolyLine> smoothLine = vtkSmartPointer<vtkPolyLine>::New();
1809  smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints());
1810 
1811  for (int j=0; j<smoothLine->GetNumberOfPoints(); j++)
1812  {
1813  smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt);
1814  vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j));
1815  }
1816  vtkSmoothCells->InsertNextCell(smoothLine);
1817  pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints();
1818  }
1819 
1820  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1821  m_FiberPolyData->SetPoints(vtkSmoothPoints);
1822  m_FiberPolyData->SetLines(vtkSmoothCells);
1823  this->SetFiberPolyData(m_FiberPolyData, true);
1824  m_FiberSampling = 10/pointDistance;
1825 }
1826 
1827 void mitk::FiberBundle::ResampleSpline(float pointDistance)
1828 {
1829  ResampleSpline(pointDistance, 0, 0, 0 );
1830 }
1831 
1832 unsigned long mitk::FiberBundle::GetNumberOfPoints()
1833 {
1834  unsigned long points = 0;
1835  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
1836  {
1837  vtkCell* cell = m_FiberPolyData->GetCell(i);
1838  points += cell->GetNumberOfPoints();
1839  }
1840  return points;
1841 }
1842 
1843 void mitk::FiberBundle::Compress(float error)
1844 {
1845  vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
1846  vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
1847 
1848  MITK_INFO << "Compressing fibers";
1849  unsigned long numRemovedPoints = 0;
1850  boost::progress_display disp(m_FiberPolyData->GetNumberOfCells());
1851 
1852  for (int i=0; i<m_FiberPolyData->GetNumberOfCells(); i++)
1853  {
1854  ++disp;
1855  vtkCell* cell = m_FiberPolyData->GetCell(i);
1856  int numPoints = cell->GetNumberOfPoints();
1857  vtkPoints* points = cell->GetPoints();
1858 
1859  // calculate curvatures
1860  std::vector< int > removedPoints; removedPoints.resize(numPoints, 0);
1861  removedPoints[0]=-1; removedPoints[numPoints-1]=-1;
1862 
1863  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
1864 
1865  bool pointFound = true;
1866  while (pointFound)
1867  {
1868  pointFound = false;
1869  double minError = error;
1870  int removeIndex = -1;
1871 
1872  for (int j=0; j<numPoints; j++)
1873  {
1874  if (removedPoints[j]==0)
1875  {
1876  double cand[3];
1877  points->GetPoint(j, cand);
1878  vnl_vector_fixed< double, 3 > candV;
1879  candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2];
1880 
1881  int validP = -1;
1882  vnl_vector_fixed< double, 3 > pred;
1883  for (int k=j-1; k>=0; k--)
1884  if (removedPoints[k]<=0)
1885  {
1886  double ref[3];
1887  points->GetPoint(k, ref);
1888  pred[0]=ref[0]; pred[1]=ref[1]; pred[2]=ref[2];
1889  validP = k;
1890  break;
1891  }
1892  int validS = -1;
1893  vnl_vector_fixed< double, 3 > succ;
1894  for (int k=j+1; k<numPoints; k++)
1895  if (removedPoints[k]<=0)
1896  {
1897  double ref[3];
1898  points->GetPoint(k, ref);
1899  succ[0]=ref[0]; succ[1]=ref[1]; succ[2]=ref[2];
1900  validS = k;
1901  break;
1902  }
1903 
1904  if (validP>=0 && validS>=0)
1905  {
1906  double a = (candV-pred).magnitude();
1907  double b = (candV-succ).magnitude();
1908  double c = (pred-succ).magnitude();
1909  double s=0.5*(a+b+c);
1910  double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c)));
1911 
1912  if (hc<minError)
1913  {
1914  removeIndex = j;
1915  minError = hc;
1916  pointFound = true;
1917  }
1918  }
1919  }
1920  }
1921 
1922  if (pointFound)
1923  {
1924  removedPoints[removeIndex] = 1;
1925  numRemovedPoints++;
1926  }
1927  }
1928 
1929  for (int j=0; j<numPoints; j++)
1930  {
1931  if (removedPoints[j]<=0)
1932  {
1933  double cand[3];
1934  points->GetPoint(j, cand);
1935  vtkIdType id = vtkNewPoints->InsertNextPoint(cand);
1936  container->GetPointIds()->InsertNextId(id);
1937  }
1938  }
1939 
1940  vtkNewCells->InsertNextCell(container);
1941  }
1942 
1943  if (vtkNewCells->GetNumberOfCells()>0)
1944  {
1945  MITK_INFO << "Removed points: " << numRemovedPoints;
1946  m_FiberPolyData = vtkSmartPointer<vtkPolyData>::New();
1947  m_FiberPolyData->SetPoints(vtkNewPoints);
1948  m_FiberPolyData->SetLines(vtkNewCells);
1949  this->SetFiberPolyData(m_FiberPolyData, true);
1950  }
1951 }
1952 
1953 // reapply selected colorcoding in case polydata structure has changed
1954 bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps)
1955 {
1956  if (fib==nullptr)
1957  {
1958  MITK_INFO << "Reference bundle is NULL!";
1959  return false;
1960  }
1961 
1962  if (m_NumFibers!=fib->GetNumFibers())
1963  {
1964  MITK_INFO << "Unequal number of fibers!";
1965  MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers();
1966  return false;
1967  }
1968 
1969  for (int i=0; i<m_NumFibers; i++)
1970  {
1971  vtkCell* cell = m_FiberPolyData->GetCell(i);
1972  int numPoints = cell->GetNumberOfPoints();
1973  vtkPoints* points = cell->GetPoints();
1974 
1975  vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i);
1976  int numPoints2 = cell2->GetNumberOfPoints();
1977  vtkPoints* points2 = cell2->GetPoints();
1978 
1979  if (numPoints2!=numPoints)
1980  {
1981  MITK_INFO << "Unequal number of points in fiber " << i << "!";
1982  MITK_INFO << numPoints2 << " vs. " << numPoints;
1983  return false;
1984  }
1985 
1986  for (int j=0; j<numPoints; j++)
1987  {
1988  double* p1 = points->GetPoint(j);
1989  double* p2 = points2->GetPoint(j);
1990  if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps)
1991  {
1992  MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!";
1993  MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2];
1994  MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2];
1995  return false;
1996  }
1997  }
1998  }
1999 
2000  return true;
2001 }
2002 
2003 /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */
2004 void mitk::FiberBundle::UpdateOutputInformation()
2005 {
2006 
2007 }
2008 void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion()
2009 {
2010 
2011 }
2012 bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion()
2013 {
2014  return false;
2015 }
2016 bool mitk::FiberBundle::VerifyRequestedRegion()
2017 {
2018  return true;
2019 }
2020 void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* )
2021 {
2022 
2023 }
mitk::FiberBundle::SubtractBundle
FiberBundle::Pointer SubtractBundle(FiberBundle *fib)
Definition: mitkFiberBundle.cpp:189
mitk::FiberBundle::RemoveLongFibers
bool RemoveLongFibers(float lengthInMM)
Definition: mitkFiberBundle.cpp:1715
mitk::FiberBundle::ResampleSpline
void ResampleSpline(float pointDistance=1)
Definition: mitkFiberBundle.cpp:1827
mitk::DataStorage
Data management class that handles 'was created by' relations.
Definition: mitkDataStorage.h:47
mitk::ImagePixelReadAccessor
Gives locked and index-based read access for a particular image part. The class provides several set-...
Definition: mitkImagePixelReadAccessor.h:36
mitk::FiberBundle::FIBER_ID_ARRAY
static const char * FIBER_ID_ARRAY
Definition: mitkFiberBundle.h:51
mitk::FiberBundle::ScaleFibers
void ScaleFibers(double x, double y, double z, bool subtractCenter=true)
Definition: mitkFiberBundle.cpp:1410
mitk::FiberBundle::SetFiberOpacity
void SetFiberOpacity(vtkDoubleArray *FAValArray)
Definition: mitkFiberBundle.cpp:537
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
mitk::FiberBundle::SetRequestedRegion
virtual void SetRequestedRegion(const itk::DataObject *) override
Set the requested region from this data object to match the requested region of the data object passe...
Definition: mitkFiberBundle.cpp:2020
mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion
virtual void SetRequestedRegionToLargestPossibleRegion() override
Set the RequestedRegion to the LargestPossibleRegion.
Definition: mitkFiberBundle.cpp:2008
mitk::FiberBundle::TransformMatrix
itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz)
Definition: mitkFiberBundle.cpp:1217
mitk::FiberBundle::UpdateOutputInformation
virtual void UpdateOutputInformation() override
Definition: mitkFiberBundle.cpp:2004
mitk::DataStorage::GetDerivations
virtual SetOfObjects::ConstPointer GetDerivations(const mitk::DataNode *node, const NodePredicateBase *condition=nullptr, bool onlyDirectDerivations=true) const =0
returns a set of derived objects for a given node.
mitk::FiberBundle::SetFiberWeights
void SetFiberWeights(float newWeight)
Definition: mitkFiberBundle.cpp:1174
mitk::FiberBundle::RemoveDir
void RemoveDir(vnl_vector_fixed< double, 3 > dir, double threshold)
Definition: mitkFiberBundle.cpp:1520
MITK_DEBUG
#define MITK_DEBUG
Definition: mitkLogMacros.h:26
std
STL namespace.
mitk::FiberBundle::GetItkPoint
itk::Point< float, 3 > GetItkPoint(double point[3])
Definition: mitkFiberBundle.cpp:260
mitk::FiberBundle::RemoveShortFibers
bool RemoveShortFibers(float lengthInMM)
Definition: mitkFiberBundle.cpp:1663
mitk::DataNode::GetData
BaseData * GetData() const
Get the data object (instance of BaseData, e.g., an Image) managed by this DataNode.
Definition: mitkDataNode.cpp:48
mitkPixelTypeMultiplex2
#define mitkPixelTypeMultiplex2(function, ptype, param1, param2)
Definition: mitkPixelTypeMultiplex.h:69
mitk::Geometry3D::New
static Pointer New()
mitk::FiberBundle::GetFiberWeight
float GetFiberWeight(unsigned int fiber)
Definition: mitkFiberBundle.cpp:1169
mitk::FiberBundle::SetFiberWeight
void SetFiberWeight(unsigned int fiber, float weight)
Definition: mitkFiberBundle.cpp:1194
mitk::FiberBundle::ItkUcharImgType
itk::Image< unsigned char, 3 > ItkUcharImgType
Definition: mitkFiberBundle.h:48
mitk::FiberBundle::TransformFibers
void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz)
Definition: mitkFiberBundle.cpp:1288
mitk::FiberBundle::RotateAroundAxis
void RotateAroundAxis(double x, double y, double z)
Definition: mitkFiberBundle.cpp:1350
mitk::ConstPointer
itk::SmartPointer< const Self > ConstPointer
Definition: mitkRenderingManager.h:389
mitk::FiberBundle::AddBundle
FiberBundle::Pointer AddBundle(FiberBundle *fib)
Definition: mitkFiberBundle.cpp:119
mitk::FiberBundle::FiberBundle
FiberBundle(vtkPolyData *fiberPolyData=nullptr)
Definition: mitkFiberBundle.cpp:53
mitk::FiberBundle::ApplyCurvatureThreshold
bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers)
Definition: mitkFiberBundle.cpp:1583
mitk::FiberBundle::GetDeepCopy
mitk::FiberBundle::Pointer GetDeepCopy()
Definition: mitkFiberBundle.cpp:76
MITK_WARN
#define MITK_WARN
Definition: mitkLogMacros.h:23
mitk::FiberBundle::TransformPoint
itk::Point< float, 3 > TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz)
Definition: mitkFiberBundle.cpp:1248
mitk::LookupTable::New
static Pointer New()
mitk::FiberBundle::GetFiberPolyData
vtkSmartPointer< vtkPolyData > GetFiberPolyData() const
Definition: mitkFiberBundle.cpp:292
mitk::FiberBundle
Base Class for Fiber Bundles;.
Definition: mitkFiberBundle.h:44
max
static T max(T x, T y)
Definition: svm.cpp:70
mitk::BaseGeometry::GetCenter
Point3D GetCenter() const
Get the center of the bounding-box in mm.
Definition: mitkBaseGeometry.cpp:268
mitk::FiberBundle::VerifyRequestedRegion
virtual bool VerifyRequestedRegion() override
Verify that the RequestedRegion is within the LargestPossibleRegion.
Definition: mitkFiberBundle.cpp:2016
mitk::FiberBundle::ExtractFiberIdSubset
std::vector< long > ExtractFiberIdSubset(DataNode *roi, DataStorage *storage)
Definition: mitkFiberBundle.cpp:906
mitk::FiberBundle::GetNumberOfPoints
unsigned long GetNumberOfPoints()
Definition: mitkFiberBundle.cpp:1832
mitk::FiberBundle::TranslateFibers
void TranslateFibers(double x, double y, double z)
Definition: mitkFiberBundle.cpp:1455
min
static T min(T x, T y)
Definition: svm.cpp:67
mitk::FiberBundle::ExtractFiberSubset
FiberBundle::Pointer ExtractFiberSubset(DataNode *roi, DataStorage *storage)
Definition: mitkFiberBundle.cpp:893
mitk::PlanarFigure
Base-class for geometric planar (2D) figures, such as lines, circles, rectangles, polygons...
Definition: mitkPlanarFigure.h:54
mitk::FiberBundle::ColorFibersByScalarMap
void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity)
Definition: mitkFiberBundle.cpp:557
mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion
virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override
Determine whether the RequestedRegion is outside of the BufferedRegion.
Definition: mitkFiberBundle.cpp:2012
mitk::FiberBundle::SetFiberColors
void SetFiberColors(vtkSmartPointer< vtkUnsignedCharArray > fiberColors)
Definition: mitkFiberBundle.cpp:1199
mitk::eps
MITKCORE_EXPORT const ScalarType eps
Definition: mitkNumericConstants.cpp:21
mitk::FiberBundle::MirrorFibers
void MirrorFibers(unsigned int axis)
Definition: mitkFiberBundle.cpp:1485
mitk::FiberBundle::RemoveFibersOutside
FiberBundle::Pointer RemoveFibersOutside(ItkUcharImgType *mask, bool invert=false)
Definition: mitkFiberBundle.cpp:815
mitk::FiberBundle::Equals
bool Equals(FiberBundle *fib, double eps=0.0001)
Definition: mitkFiberBundle.cpp:1954
mitk::FiberBundle::Compress
void Compress(float error=0.0)
Definition: mitkFiberBundle.cpp:1843
mitk::FiberBundle::GetNumFibers
virtual int GetNumFibers()
mitk::FiberBundle::SetFiberPolyData
void SetFiberPolyData(vtkSmartPointer< vtkPolyData >, bool updateGeometry=true)
Definition: mitkFiberBundle.cpp:272
mitk::BaseGeometry
BaseGeometry Describes the geometry of a data object.
Definition: mitkBaseGeometry.h:101
mitk::DataNode
Class for nodes of the DataTree.
Definition: mitkDataNode.h:66
mitk::PixelType
Class for defining the data type of pixels.
Definition: mitkPixelType.h:55
mitk::FiberBundle::GeneratePolyDataByIds
vtkSmartPointer< vtkPolyData > GeneratePolyDataByIds(std::vector< long >)
Definition: mitkFiberBundle.cpp:84
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.