Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
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.