1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
18 #ifndef __mitkOdfVtkMapper2D_txx__
19 #define __mitkOdfVtkMapper2D_txx__
21 #include "mitkOdfVtkMapper2D.h"
22 #include "mitkDataNode.h"
23 #include "mitkBaseRenderer.h"
24 #include "mitkMatrixConvert.h"
25 #include "mitkGeometry3D.h"
26 #include "mitkTimeGeometry.h"
27 #include "mitkOdfNormalizationMethodProperty.h"
28 #include "mitkOdfScaleByProperty.h"
29 #include "mitkProperties.h"
30 #include "mitkTensorImage.h"
32 #include "vtkSphereSource.h"
33 #include "vtkPropCollection.h"
34 #include "vtkMaskedGlyph3D.h"
35 #include "vtkGlyph2D.h"
36 #include "vtkGlyph3D.h"
37 #include "vtkMaskedProgrammableGlyphFilter.h"
38 #include "vtkImageData.h"
39 #include "vtkLinearTransform.h"
40 #include "vtkCamera.h"
41 #include "vtkPointData.h"
42 #include "vtkTransformPolyDataFilter.h"
43 #include "vtkTransform.h"
44 #include "vtkOdfSource.h"
45 #include "vtkDoubleArray.h"
46 #include "vtkLookupTable.h"
47 #include "vtkProperty.h"
48 #include "vtkPolyDataNormals.h"
50 #include "vtkLightCollection.h"
52 #include "vtkFloatArray.h"
53 #include "vtkDelaunay2D.h"
54 #include "vtkMapper.h"
55 #include <vtkInformationVector.h>
56 #include <vtkInformation.h>
57 #include "vtkRenderer.h"
59 #include "itkOrientationDistributionFunction.h"
61 #include "itkFixedArray.h"
64 #include "vtkOpenGLRenderer.h"
66 #define _USE_MATH_DEFINES
72 template<class T, int N>
73 vtkSmartPointer<vtkTransform> mitk::OdfVtkMapper2D<T,N>::m_OdfTransform = vtkSmartPointer<vtkTransform>::New();
75 template<class T, int N>
76 vtkSmartPointer<vtkOdfSource> mitk::OdfVtkMapper2D<T,N>::m_OdfSource = vtkSmartPointer<vtkOdfSource>::New();
78 template<class T, int N>
79 float mitk::OdfVtkMapper2D<T,N>::m_Scaling;
81 template<class T, int N>
82 int mitk::OdfVtkMapper2D<T,N>::m_Normalization;
84 template<class T, int N>
85 int mitk::OdfVtkMapper2D<T,N>::m_ScaleBy;
87 template<class T, int N>
88 float mitk::OdfVtkMapper2D<T,N>::m_IndexParam1;
90 template<class T, int N>
91 float mitk::OdfVtkMapper2D<T,N>::m_IndexParam2;
93 template<class T, int N>
94 bool mitk::OdfVtkMapper2D<T, N>::m_toggleTensorEllipsoidView = false;
96 template<class T, int N>
97 bool mitk::OdfVtkMapper2D<T, N>::m_toggleColourisationMode = false;
99 template<class T, int N>
100 bool mitk::OdfVtkMapper2D<T, N>::m_toggleGlyphPlacementMode = true;
102 template<class T, int N>
103 vtkSmartPointer<vtkDoubleArray> mitk::OdfVtkMapper2D<T, N>::m_colourScalars = nullptr;
105 #define ODF_MAPPER_PI M_PI
108 template<class T, int N>
109 mitk::OdfVtkMapper2D<T,N>::LocalStorage::LocalStorage()
111 m_PropAssemblies.push_back(vtkPropAssembly::New());
112 m_PropAssemblies.push_back(vtkPropAssembly::New());
113 m_PropAssemblies.push_back(vtkPropAssembly::New());
115 m_OdfsPlanes.push_back(vtkAppendPolyData::New());
116 m_OdfsPlanes.push_back(vtkAppendPolyData::New());
117 m_OdfsPlanes.push_back(vtkAppendPolyData::New());
119 m_OdfsPlanes[0]->AddInputData(vtkPolyData::New());
120 m_OdfsPlanes[1]->AddInputData(vtkPolyData::New());
121 m_OdfsPlanes[2]->AddInputData(vtkPolyData::New());
123 m_OdfsActors.push_back(vtkActor::New());
124 m_OdfsActors.push_back(vtkActor::New());
125 m_OdfsActors.push_back(vtkActor::New());
127 m_OdfsActors[0]->GetProperty()->SetInterpolationToGouraud();
128 m_OdfsActors[1]->GetProperty()->SetInterpolationToGouraud();
129 m_OdfsActors[2]->GetProperty()->SetInterpolationToGouraud();
131 m_OdfsMappers.push_back(vtkPolyDataMapper::New());
132 m_OdfsMappers.push_back(vtkPolyDataMapper::New());
133 m_OdfsMappers.push_back(vtkPolyDataMapper::New());
135 vtkLookupTable *lut = vtkLookupTable::New();
137 m_OdfsMappers[0]->SetLookupTable(lut);
138 m_OdfsMappers[1]->SetLookupTable(lut);
139 m_OdfsMappers[2]->SetLookupTable(lut);
141 m_OdfsActors[0]->SetMapper(m_OdfsMappers[0]);
142 m_OdfsActors[1]->SetMapper(m_OdfsMappers[1]);
143 m_OdfsActors[2]->SetMapper(m_OdfsMappers[2]);
146 template<class T, int N>
147 mitk::OdfVtkMapper2D<T,N>
151 m_Planes.push_back(vtkPlane::New());
152 m_Planes.push_back(vtkPlane::New());
153 m_Planes.push_back(vtkPlane::New());
155 m_Cutters.push_back(vtkCutter::New());
156 m_Cutters.push_back(vtkCutter::New());
157 m_Cutters.push_back(vtkCutter::New());
159 m_Cutters[0]->SetCutFunction( m_Planes[0] );
160 m_Cutters[0]->GenerateValues( 1, 0, 1 );
162 m_Cutters[1]->SetCutFunction( m_Planes[1] );
163 m_Cutters[1]->GenerateValues( 1, 0, 1 );
165 m_Cutters[2]->SetCutFunction( m_Planes[2] );
166 m_Cutters[2]->GenerateValues( 1, 0, 1 );
168 // Windowing the cutted planes in direction 1
169 m_ThickPlanes1.push_back(vtkThickPlane::New());
170 m_ThickPlanes1.push_back(vtkThickPlane::New());
171 m_ThickPlanes1.push_back(vtkThickPlane::New());
173 m_Clippers1.push_back(vtkClipPolyData::New());
174 m_Clippers1.push_back(vtkClipPolyData::New());
175 m_Clippers1.push_back(vtkClipPolyData::New());
177 m_Clippers1[0]->SetClipFunction( m_ThickPlanes1[0] );
178 m_Clippers1[1]->SetClipFunction( m_ThickPlanes1[1] );
179 m_Clippers1[2]->SetClipFunction( m_ThickPlanes1[2] );
181 // Windowing the cutted planes in direction 2
182 m_ThickPlanes2.push_back(vtkThickPlane::New());
183 m_ThickPlanes2.push_back(vtkThickPlane::New());
184 m_ThickPlanes2.push_back(vtkThickPlane::New());
186 m_Clippers2.push_back(vtkClipPolyData::New());
187 m_Clippers2.push_back(vtkClipPolyData::New());
188 m_Clippers2.push_back(vtkClipPolyData::New());
190 m_Clippers2[0]->SetClipFunction( m_ThickPlanes2[0] );
191 m_Clippers2[1]->SetClipFunction( m_ThickPlanes2[1] );
192 m_Clippers2[2]->SetClipFunction( m_ThickPlanes2[2] );
194 m_ShowMaxNumber = 500;
197 template<class T, int N>
198 mitk::OdfVtkMapper2D<T,N>
204 template<class T, int N>
205 mitk::Image* mitk::OdfVtkMapper2D<T,N>
208 return static_cast<mitk::Image * > ( m_DataNode->GetData() );
211 template<class T, int N>
212 vtkProp* mitk::OdfVtkMapper2D<T,N>
213 ::GetVtkProp(mitk::BaseRenderer* renderer)
215 LocalStorage *localStorage = m_LSH.GetLocalStorage(renderer);
216 return localStorage->m_PropAssemblies[GetIndex(renderer)];
219 template<class T, int N>
220 int mitk::OdfVtkMapper2D<T,N>
221 ::GetIndex(mitk::BaseRenderer* renderer)
223 if(!strcmp(renderer->GetName(),"stdmulti.widget1"))
226 if(!strcmp(renderer->GetName(),"stdmulti.widget2"))
229 if(!strcmp(renderer->GetName(),"stdmulti.widget3"))
235 template<class T, int N>
236 void mitk::OdfVtkMapper2D<T,N>
237 ::GlyphMethod(void *arg)
239 vtkMaskedProgrammableGlyphFilter* pfilter=(vtkMaskedProgrammableGlyphFilter*)arg;
242 double debugpoint[3];
243 pfilter->GetPoint(point);
244 pfilter->GetPoint(debugpoint);
246 itk::Point<double,3> p(point);
247 Vector3D spacing = pfilter->GetGeometry()->GetSpacing();
253 pfilter->GetGeometry()->IndexToWorld( p, p2 );
258 vtkPointData* data = pfilter->GetPointData();
259 vtkDataArray* odfvals = data->GetArray("vector");
260 vtkIdType id = pfilter->GetPointId();
261 m_OdfTransform->Identity();
262 m_OdfTransform->Translate(point[0],point[1],point[2]);
264 typedef itk::OrientationDistributionFunction<float,N> OdfType;
267 if( odfvals->GetNumberOfComponents()==6 and not m_toggleTensorEllipsoidView )
269 float tensorelems[6] = {
270 (float)odfvals->GetComponent(id,0),
271 (float)odfvals->GetComponent(id,1),
272 (float)odfvals->GetComponent(id,2),
273 (float)odfvals->GetComponent(id,3),
274 (float)odfvals->GetComponent(id,4),
275 (float)odfvals->GetComponent(id,5),
277 itk::DiffusionTensor3D<float> tensor(tensorelems);
278 odf.InitFromTensor(tensor);
280 if( (pfilter-> GetColorMode()) == 0 ) // m_ColourisationMode; VTK_COLOR_BY_INPUT=0, VTK_COLOR_BY_SOURCE=1.
281 { /// \brief Colourisation of glyph like in MitkWorkbench's dti visualisation, r,g,b,a=x,y,z,fa of main direction of diffusion.
282 // unsigned char rgba[4] = {0,0,0,0};
283 double rgba[4] = {0,0,0,0};
284 rgba[3] = fabs(tensor.GetFractionalAnisotropy()); //* 255); // fa = FractionalAnisotropy => Helligkeit = Value, and alpha.
285 typename itk::DiffusionTensor3D<T>::EigenValuesArrayType eigenValues;
286 typename itk::DiffusionTensor3D<T>::EigenVectorsMatrixType eigenVectors;
287 tensor.ComputeEigenAnalysis( eigenValues, eigenVectors ); // normalized eigenvectors as rows in ascending order.
288 rgba[0] = (fabs( eigenVectors(2, 0) ) * rgba[3]);
289 rgba[1] = (fabs( eigenVectors(2, 1) ) * rgba[3]);
290 rgba[2] = (fabs( eigenVectors(2, 2) ) * rgba[3]);
291 //#define __IMG_DAT_ITEM__CLIP_00_FF__(val) (val) = ( (val) < 0 ) ? ( 0 ) : ( ( (val) > 255 ) ? ( 255 ) : ( (val) ) );
292 #define __IMG_DAT_ITEM__CLIP_0_1__(val) (val) = ( (val) < 0.0 ) ? ( 0.0 ) : ( ( (val) > 1.0 ) ? ( 1.0 ) : ( (val) ) );
293 __IMG_DAT_ITEM__CLIP_0_1__(rgba[0]);
294 __IMG_DAT_ITEM__CLIP_0_1__(rgba[1]);
295 __IMG_DAT_ITEM__CLIP_0_1__(rgba[2]);
296 __IMG_DAT_ITEM__CLIP_0_1__(rgba[3]);
297 m_colourScalars-> InsertNextTuple4( rgba[0], rgba[1], rgba[2], rgba[3] ); // Allocates space automaticaly.
299 data-> AddArray( m_colourScalars );
300 data-> SetActiveScalars( "GLYPH_COLORS" );
304 data-> SetActiveScalars( "vector" );
305 data-> RemoveArray( "GLYPH_COLORS" );
309 else if( odfvals->GetNumberOfComponents()==6 and m_toggleTensorEllipsoidView )
311 float tensorElements[6] = {
312 (float)odfvals->GetComponent(id, 0), (float)odfvals->GetComponent(id, 1), (float)odfvals->GetComponent(id, 2),
313 (float)odfvals->GetComponent(id, 3), (float)odfvals->GetComponent(id, 4), (float)odfvals->GetComponent(id, 5),
315 itk::DiffusionTensor3D<float> tensor( tensorElements );
316 odf.InitFromEllipsoid( tensor );
318 if( (pfilter-> GetColorMode()) == 0 ) // m_ColourisationMode; VTK_COLOR_BY_INPUT=0, VTK_COLOR_BY_SOURCE=1.
319 { /// \brief Colourisation of glyph like in MitkWorkbench's dti visualisation, r,g,b,a=x,y,z,fa of main direction of diffusion.
320 double rgba[4] = {0,0,0,0};
321 rgba[3] = fabs( tensor.GetFractionalAnisotropy()); // * 255 ); // fa = FractionalAnisotropy => Helligkeit = Value
322 typename itk::DiffusionTensor3D<T>::EigenValuesArrayType eigenValues;
323 typename itk::DiffusionTensor3D<T>::EigenVectorsMatrixType eigenVectors;
324 tensor.ComputeEigenAnalysis( eigenValues, eigenVectors ); // normalized eigenvectors as rows in ascending order.
325 rgba[0] = (fabs( eigenVectors(2, 0) ) * rgba[3]);
326 rgba[1] = (fabs( eigenVectors(2, 1) ) * rgba[3]);
327 rgba[2] = (fabs( eigenVectors(2, 2) ) * rgba[3]);
328 __IMG_DAT_ITEM__CLIP_0_1__(rgba[0]);
329 __IMG_DAT_ITEM__CLIP_0_1__(rgba[1]);
330 __IMG_DAT_ITEM__CLIP_0_1__(rgba[2]);
331 __IMG_DAT_ITEM__CLIP_0_1__(rgba[3]);
333 m_colourScalars-> InsertNextTuple4( rgba[0], rgba[1], rgba[2], rgba[3] );
334 data-> CopyScalarsOn();
335 data-> SetCopyAttribute(id, 1);
336 data-> SetCopyScalars(id);
338 data-> SetDebug('1');
341 pfilter-> GetOutput()-> DebugOn();
342 pfilter-> GetOutput()-> GetPointData()-> CopyScalarsOn();
344 data-> AddArray( m_colourScalars );
345 data-> SetActiveScalars( "GLYPH_COLORS" );
350 std::cout << "<data>\n";
351 data-> PrintSelf(std::cout, indent.GetNextIndent());
352 std::cout << "</data>\n<pfilter>\n";
353 pfilter-> PrintSelf(std::cout, indent.GetNextIndent());
354 std::cout << pfilter->GetOutput()->GetClassName() << "\n";
355 pfilter-> GetOutput()-> PrintSelf(std::cout, indent.GetNextIndent());
356 pfilter-> GetOutput()->GetPointData()-> PrintSelf(std::cout, indent.GetNextIndent());
357 pfilter-> GetOutput()->GetPointData()->GetScalars()-> PrintSelf(std::cout, indent.GetNextIndent());
358 std::cout << "</pfilter>";
365 data-> SetActiveScalars( "vector" );
366 data-> RemoveArray( "GLYPH_COLORS" );
372 for(int i=0; i<N; i++)
373 odf[i] = (double)odfvals->GetComponent(id,i);
379 m_OdfSource->SetScale(m_Scaling);
382 m_OdfSource->SetScale(m_Scaling*odf.GetGeneralizedGFA(m_IndexParam1, m_IndexParam2));
385 m_OdfSource->SetScale(m_Scaling*odf.GetPrincipleCurvature(m_IndexParam1, m_IndexParam2, 0));
389 m_OdfSource->SetNormalization(m_Normalization);
390 m_OdfSource->SetOdf(odf);
391 m_OdfSource->Modified();
394 template<class T, int N>
395 typename mitk::OdfVtkMapper2D<T,N>::OdfDisplayGeometry mitk::OdfVtkMapper2D<T,N>
396 ::MeasureDisplayedGeometry(mitk::BaseRenderer* renderer)
398 PlaneGeometry::ConstPointer worldPlaneGeometry = renderer->GetCurrentWorldPlaneGeometry();
400 // set up the cutter orientation according to the current geometry of
401 // the renderers plane
402 double vp[ 3 ], vnormal[ 3 ];
403 Point3D point = worldPlaneGeometry->GetOrigin();
404 Vector3D normal = worldPlaneGeometry->GetNormal(); normal.Normalize();
405 vnl2vtk( point.GetVnlVector(), vp );
406 vnl2vtk( normal.GetVnlVector(), vnormal );
408 Point2D dispSizeInMM = renderer->GetViewportSizeInMM();
410 Point2D displayGeometryOriginInMM = renderer->GetOriginInMM();
412 mitk::Vector2D size = dispSizeInMM.GetVectorFromOrigin();
413 mitk::Vector2D origin = displayGeometryOriginInMM.GetVectorFromOrigin();
427 M[0] = origin[0] + size[0]/2;
428 M[1] = origin[1] + size[1]/2;
431 L[1] = origin[1] + size[1]/2;
433 O[0] = origin[0] + size[0]/2;
434 O[1] = origin[1] + size[1];
436 mitk::Point2D point1;
437 point1[0] = M[0]; point1[1] = M[1];
439 renderer->GetCurrentWorldPlaneGeometry()->Map(point1, M3D);
441 point1[0] = L[0]; point1[1] = L[1];
443 renderer->GetCurrentWorldPlaneGeometry()->Map(point1, L3D);
445 point1[0] = O[0]; point1[1] = O[1];
447 renderer->GetCurrentWorldPlaneGeometry()->Map(point1, O3D);
449 double d1 = sqrt((M3D[0]-L3D[0])*(M3D[0]-L3D[0])
450 + (M3D[1]-L3D[1])*(M3D[1]-L3D[1])
451 + (M3D[2]-L3D[2])*(M3D[2]-L3D[2]));
452 double d2 = sqrt((M3D[0]-O3D[0])*(M3D[0]-O3D[0])
453 + (M3D[1]-O3D[1])*(M3D[1]-O3D[1])
454 + (M3D[2]-O3D[2])*(M3D[2]-O3D[2]));
455 double d = d1>d2 ? d1 : d2;
458 OdfDisplayGeometry retval;
459 retval.vp[0] = vp[0];
460 retval.vp[1] = vp[1];
461 retval.vp[2] = vp[2];
462 retval.vnormal[0] = vnormal[0];
463 retval.vnormal[1] = vnormal[1];
464 retval.vnormal[2] = vnormal[2];
465 retval.normal[0] = normal[0];
466 retval.normal[1] = normal[1];
467 retval.normal[2] = normal[2];
471 retval.M3D[0] = M3D[0];
472 retval.M3D[1] = M3D[1];
473 retval.M3D[2] = M3D[2];
474 retval.L3D[0] = L3D[0];
475 retval.L3D[1] = L3D[1];
476 retval.L3D[2] = L3D[2];
477 retval.O3D[0] = O3D[0];
478 retval.O3D[1] = O3D[1];
479 retval.O3D[2] = O3D[2];
481 retval.vp_original[0] = vp[0];
482 retval.vp_original[1] = vp[1];
483 retval.vp_original[2] = vp[2];
484 retval.vnormal_original[0] = vnormal[0];
485 retval.vnormal_original[1] = vnormal[1];
486 retval.vnormal_original[2] = vnormal[2];
487 retval.size[0] = size[0];
488 retval.size[1] = size[1];
489 retval.origin[0] = origin[0];
490 retval.origin[1] = origin[1];
495 template<class T, int N>
496 void mitk::OdfVtkMapper2D<T,N>
497 ::Slice(mitk::BaseRenderer* renderer, OdfDisplayGeometry dispGeo)
499 LocalStorage *localStorage = m_LSH.GetLocalStorage(renderer);
501 vtkLinearTransform * vtktransform =
502 this->GetDataNode()->GetVtkTransform(this->GetTimestep());
504 int index = GetIndex(renderer);
506 vtkSmartPointer<vtkTransform> inversetransform = vtkSmartPointer<vtkTransform>::New();
507 inversetransform->Identity();
508 inversetransform->Concatenate(vtktransform->GetLinearInverse());
510 ((vtkTransform*)vtktransform)->GetScale(myscale);
511 inversetransform->PostMultiply();
512 inversetransform->Scale(1*myscale[0],1*myscale[1],1*myscale[2]);
513 inversetransform->TransformPoint( dispGeo.vp, dispGeo.vp );
514 inversetransform->TransformNormalAtPoint( dispGeo.vp, dispGeo.vnormal, dispGeo.vnormal );
516 // vtk works in axis align coords
517 // thus the normal also must be axis align, since
518 // we do not allow arbitrary cutting through volume
520 // vnormal should already be axis align, but in order
521 // to get rid of precision effects, we set the two smaller
522 // components to zero here
524 m_VtkImage->GetDimensions(dims);
526 m_VtkImage->GetSpacing(spac);
527 if(fabs(dispGeo.vnormal[0]) > fabs(dispGeo.vnormal[1])
528 && fabs(dispGeo.vnormal[0]) > fabs(dispGeo.vnormal[2]) )
530 if(fabs(dispGeo.vp[0]/spac[0]) < 0.4)
531 dispGeo.vp[0] = 0.4*spac[0];
532 if(fabs(dispGeo.vp[0]/spac[0]) > (dims[0]-1)-0.4)
533 dispGeo.vp[0] = ((dims[0]-1)-0.4)*spac[0];
534 dispGeo.vnormal[1] = 0;
535 dispGeo.vnormal[2] = 0;
538 if(fabs(dispGeo.vnormal[1]) > fabs(dispGeo.vnormal[0]) && fabs(dispGeo.vnormal[1]) > fabs(dispGeo.vnormal[2]) )
540 if(fabs(dispGeo.vp[1]/spac[1]) < 0.4)
541 dispGeo.vp[1] = 0.4*spac[1];
542 if(fabs(dispGeo.vp[1]/spac[1]) > (dims[1]-1)-0.4)
543 dispGeo.vp[1] = ((dims[1]-1)-0.4)*spac[1];
544 dispGeo.vnormal[0] = 0;
545 dispGeo.vnormal[2] = 0;
548 if(fabs(dispGeo.vnormal[2]) > fabs(dispGeo.vnormal[1]) && fabs(dispGeo.vnormal[2]) > fabs(dispGeo.vnormal[0]) )
550 if(fabs(dispGeo.vp[2]/spac[2]) < 0.4)
551 dispGeo.vp[2] = 0.4*spac[2];
552 if(fabs(dispGeo.vp[2]/spac[2]) > (dims[2]-1)-0.4)
553 dispGeo.vp[2] = ((dims[2]-1)-0.4)*spac[2];
554 dispGeo.vnormal[0] = 0;
555 dispGeo.vnormal[1] = 0;
559 m_Planes[index]->SetTransform( (vtkAbstractTransform*)NULL );
560 m_Planes[index]->SetOrigin( dispGeo.vp );
561 m_Planes[index]->SetNormal( dispGeo.vnormal );
563 vtkSmartPointer<vtkPoints> points;
564 vtkSmartPointer<vtkPoints> tmppoints;
565 vtkSmartPointer<vtkPolyData> polydata;
566 vtkSmartPointer<vtkFloatArray> pointdata;
567 vtkSmartPointer<vtkDelaunay2D> delaunay;
568 vtkSmartPointer<vtkPolyData> cuttedPlane;
570 // the cutter only works if we do not have a 2D-image
571 // or if we have a 2D-image and want to see the whole image.
573 // for side views of 2D-images, we need some special treatment
574 if(!( (dims[0] == 1 && dispGeo.vnormal[0] != 0) ||
575 (dims[1] == 1 && dispGeo.vnormal[1] != 0) ||
576 (dims[2] == 1 && dispGeo.vnormal[2] != 0) ))
578 m_Cutters[index]->SetCutFunction( m_Planes[index] );
579 m_Cutters[index]->SetInputData( m_VtkImage );
580 m_Cutters[index]->Update();
581 cuttedPlane = m_Cutters[index]->GetOutput();
585 // cutting of a 2D-Volume does not work,
586 // so we have to build up our own polydata object
587 cuttedPlane = vtkPolyData::New();
588 points = vtkPoints::New();
589 points->SetNumberOfPoints(m_VtkImage->GetNumberOfPoints());
590 for(int i=0; i<m_VtkImage->GetNumberOfPoints(); i++)
592 points->SetPoint(i, m_VtkImage->GetPoint(i));
594 cuttedPlane->SetPoints(points);
596 pointdata = vtkFloatArray::New();
597 int comps = m_VtkImage->GetPointData()->GetScalars()->GetNumberOfComponents();
598 pointdata->SetNumberOfComponents(comps);
599 int tuples = m_VtkImage->GetPointData()->GetScalars()->GetNumberOfTuples();
600 pointdata->SetNumberOfTuples(tuples);
601 for(int i=0; i<tuples; i++)
603 pointdata->SetTuple( i, m_VtkImage->GetPointData()->GetScalars()->GetTuple(i) );
605 pointdata->SetName( "vector" );
606 cuttedPlane->GetPointData()->AddArray(pointdata);
611 nZero1 = 1; nZero2 = 2;
615 nZero1 = 0; nZero2 = 2;
619 nZero1 = 0; nZero2 = 1;
622 tmppoints = vtkPoints::New();
623 for(int j=0; j<m_VtkImage->GetNumberOfPoints(); j++){
625 m_VtkImage->GetPoint(j,pt);
626 tmppoints->InsertNextPoint(pt[nZero1],pt[nZero2],0);
629 polydata = vtkPolyData::New();
630 polydata->SetPoints( tmppoints );
631 delaunay = vtkDelaunay2D::New();
632 delaunay->SetInputData( polydata );
634 vtkCellArray* polys = delaunay->GetOutput()->GetPolys();
635 cuttedPlane->SetPolys(polys);
638 if(cuttedPlane->GetNumberOfPoints())
641 inversetransform = vtkTransform::New();
642 inversetransform->Identity();
643 inversetransform->Concatenate(vtktransform->GetLinearInverse());
645 ((vtkTransform*)vtktransform)->GetScale(myscale);
646 inversetransform->PostMultiply();
647 inversetransform->Scale(1*myscale[0],1*myscale[1],1*myscale[2]);
649 dispGeo.vnormal[0] = dispGeo.M3D[0]-dispGeo.O3D[0];
650 dispGeo.vnormal[1] = dispGeo.M3D[1]-dispGeo.O3D[1];
651 dispGeo.vnormal[2] = dispGeo.M3D[2]-dispGeo.O3D[2];
652 vtkMath::Normalize(dispGeo.vnormal);
653 dispGeo.vp[0] = dispGeo.M3D[0];
654 dispGeo.vp[1] = dispGeo.M3D[1];
655 dispGeo.vp[2] = dispGeo.M3D[2];
657 inversetransform->TransformPoint( dispGeo.vp, dispGeo.vp );
658 inversetransform->TransformNormalAtPoint( dispGeo.vp, dispGeo.vnormal, dispGeo.vnormal );
660 m_ThickPlanes1[index]->count = 0;
661 m_ThickPlanes1[index]->SetTransform((vtkAbstractTransform*)NULL );
662 m_ThickPlanes1[index]->SetPose( dispGeo.vnormal, dispGeo.vp );
663 m_ThickPlanes1[index]->SetThickness(dispGeo.d2);
664 m_Clippers1[index]->SetClipFunction( m_ThickPlanes1[index] );
665 m_Clippers1[index]->SetInputData( cuttedPlane );
666 m_Clippers1[index]->SetInsideOut(1);
667 m_Clippers1[index]->Update();
669 dispGeo.vnormal[0] = dispGeo.M3D[0]-dispGeo.L3D[0];
670 dispGeo.vnormal[1] = dispGeo.M3D[1]-dispGeo.L3D[1];
671 dispGeo.vnormal[2] = dispGeo.M3D[2]-dispGeo.L3D[2];
672 vtkMath::Normalize(dispGeo.vnormal);
673 dispGeo.vp[0] = dispGeo.M3D[0];
674 dispGeo.vp[1] = dispGeo.M3D[1];
675 dispGeo.vp[2] = dispGeo.M3D[2];
677 inversetransform->TransformPoint( dispGeo.vp, dispGeo.vp );
678 inversetransform->TransformNormalAtPoint( dispGeo.vp, dispGeo.vnormal, dispGeo.vnormal );
680 m_ThickPlanes2[index]->count = 0;
681 m_ThickPlanes2[index]->SetTransform((vtkAbstractTransform*)NULL );
682 m_ThickPlanes2[index]->SetPose( dispGeo.vnormal, dispGeo.vp );
683 m_ThickPlanes2[index]->SetThickness(dispGeo.d1);
684 m_Clippers2[index]->SetClipFunction( m_ThickPlanes2[index] );
685 m_Clippers2[index]->SetInputData( m_Clippers1[index]->GetOutput() );
686 m_Clippers2[index]->SetInsideOut(1);
687 m_Clippers2[index]->Update();
689 cuttedPlane = m_Clippers2[index]->GetOutput ();
691 if(cuttedPlane->GetNumberOfPoints())
693 localStorage->m_OdfsPlanes[index]->RemoveAllInputs();
695 vtkSmartPointer<vtkPolyDataNormals> normals = vtkSmartPointer<vtkPolyDataNormals>::New();
696 normals->SetInputConnection( m_OdfSource->GetOutputPort() );
697 normals->SplittingOff();
698 normals->ConsistencyOff();
699 normals->AutoOrientNormalsOff();
700 normals->ComputePointNormalsOn();
701 normals->ComputeCellNormalsOff();
702 normals->FlipNormalsOff();
703 normals->NonManifoldTraversalOff();
705 vtkSmartPointer<vtkTransformPolyDataFilter> trans = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
706 trans->SetInputConnection( normals->GetOutputPort() );
707 trans->SetTransform(m_OdfTransform);
709 vtkSmartPointer<vtkMaskedProgrammableGlyphFilter> glyphGenerator = vtkSmartPointer<vtkMaskedProgrammableGlyphFilter>::New();
710 glyphGenerator->SetMaximumNumberOfPoints(std::min(m_ShowMaxNumber,(int)cuttedPlane->GetNumberOfPoints()));
712 glyphGenerator->SetRandomMode( m_toggleGlyphPlacementMode );
714 glyphGenerator->SetUseMaskPoints(1);
715 glyphGenerator->SetSourceConnection(trans->GetOutputPort() );
716 glyphGenerator->SetInput(cuttedPlane);
718 if( not m_toggleColourisationMode )
720 glyphGenerator-> SetColorModeToColorBySource();
721 glyphGenerator-> SetInputArrayToProcess(0,0,0, vtkDataObject::FIELD_ASSOCIATION_POINTS , "vector");
722 if(m_colourScalars != nullptr) { m_colourScalars-> Initialize(); }
724 else if ( m_toggleColourisationMode )
726 m_colourScalars = vtkSmartPointer<vtkDoubleArray>::New();
727 m_colourScalars-> SetNumberOfComponents( 4 ); // red, green, blue and alpha are the 4 components per tuple.
728 m_colourScalars-> SetName( "GLYPH_COLORS" );
729 glyphGenerator-> SetColorModeToColorByInput();
730 glyphGenerator-> SetInputArrayToProcess(0,0,0, vtkDataObject::FIELD_ASSOCIATION_POINTS , "GLYPH_COLORS");
733 glyphGenerator->SetGeometry(this->GetDataNode()->GetData()->GetGeometry());
734 glyphGenerator->SetGlyphMethod(&(GlyphMethod),(void *)glyphGenerator);
738 glyphGenerator->Update();
740 catch( itk::ExceptionObject& err )
742 std::cout << err << std::endl;
745 localStorage->m_OdfsPlanes[index]->AddInputConnection(glyphGenerator->GetOutputPort());
746 localStorage->m_OdfsPlanes[index]->Update();
749 localStorage->m_PropAssemblies[index]->VisibilityOn();
750 if(localStorage->m_PropAssemblies[index]->GetParts()->IsItemPresent(localStorage->m_OdfsActors[index]))
752 localStorage->m_PropAssemblies[index]->RemovePart(localStorage->m_OdfsActors[index]);
754 localStorage->m_OdfsMappers[index]->SetInputData(localStorage->m_OdfsPlanes[index]->GetOutput());
755 localStorage->m_PropAssemblies[index]->AddPart(localStorage->m_OdfsActors[index]);
758 template<class T, int N>
759 bool mitk::OdfVtkMapper2D<T,N>
760 ::IsVisibleOdfs(mitk::BaseRenderer* renderer)
762 mitk::Image::Pointer input = const_cast<mitk::Image*>(this->GetInput());
763 const TimeGeometry *inputTimeGeometry = input->GetTimeGeometry();
764 if(inputTimeGeometry==NULL || inputTimeGeometry->CountTimeSteps()==0 || !inputTimeGeometry->IsValidTimeStep(this->GetTimestep()))
767 if(this->IsPlaneRotated(renderer))
771 switch(GetIndex(renderer))
774 GetDataNode()->GetVisibility(retval, renderer, "VisibleOdfs_T");
777 GetDataNode()->GetVisibility(retval, renderer, "VisibleOdfs_S");
780 GetDataNode()->GetVisibility(retval, renderer, "VisibleOdfs_C");
787 template<class T, int N>
788 void mitk::OdfVtkMapper2D<T,N>
789 ::MitkRenderOverlay(mitk::BaseRenderer* renderer)
791 if ( this->IsVisibleOdfs(renderer)==false )
794 if ( this->GetVtkProp(renderer)->GetVisibility() )
795 this->GetVtkProp(renderer)->RenderOverlay(renderer->GetVtkRenderer());
798 template<class T, int N>
799 void mitk::OdfVtkMapper2D<T,N>
800 ::MitkRenderOpaqueGeometry(mitk::BaseRenderer* renderer)
802 if ( this->IsVisibleOdfs( renderer )==false )
805 if ( this->GetVtkProp(renderer)->GetVisibility() )
808 this->GetVtkProp(renderer)->RenderOpaqueGeometry( renderer->GetVtkRenderer() );
812 template<class T, int N>
813 void mitk::OdfVtkMapper2D<T,N>
814 ::MitkRenderTranslucentGeometry(mitk::BaseRenderer* renderer)
816 if ( this->IsVisibleOdfs(renderer)==false )
819 if ( this->GetVtkProp(renderer)->GetVisibility() )
820 this->GetVtkProp(renderer)->RenderTranslucentPolygonalGeometry(renderer->GetVtkRenderer());
824 template<class T, int N>
825 void mitk::OdfVtkMapper2D<T,N>
826 ::Update(mitk::BaseRenderer* renderer)
829 GetDataNode()->GetVisibility(visible, renderer, "visible");
831 if ( !visible ) return;
833 mitk::Image::Pointer input = const_cast<mitk::Image*>( this->GetInput() );
834 if ( input.IsNull() ) return ;
836 std::string classname("TensorImage");
837 if(classname.compare(input->GetNameOfClass())==0)
838 m_VtkImage = dynamic_cast<mitk::TensorImage*>( this->GetInput() )->GetNonRgbVtkImageData();
840 std::string qclassname("QBallImage");
841 if(qclassname.compare(input->GetNameOfClass())==0)
842 m_VtkImage = dynamic_cast<mitk::QBallImage*>( this->GetInput() )->GetNonRgbVtkImageData();
846 // make sure, that we have point data with more than 1 component (as vectors)
847 vtkPointData* pointData = m_VtkImage->GetPointData();
848 if ( pointData == NULL )
850 itkWarningMacro( << "m_VtkImage->GetPointData() returns NULL!" );
853 if ( pointData->GetNumberOfArrays() == 0 )
855 itkWarningMacro( << "m_VtkImage->GetPointData()->GetNumberOfArrays() is 0!" );
858 else if ( pointData->GetArray(0)->GetNumberOfComponents() != N
859 && pointData->GetArray(0)->GetNumberOfComponents() != 6 /*for tensor visualization*/)
861 itkWarningMacro( << "number of components != number of directions in ODF!" );
864 else if ( pointData->GetArrayName( 0 ) == NULL )
866 m_VtkImage->GetPointData()->GetArray(0)->SetName("vector");
869 GenerateDataForRenderer(renderer);
873 itkWarningMacro( << "m_VtkImage is NULL!" );
878 template<class T, int N>
879 void mitk::OdfVtkMapper2D<T,N>
880 ::GenerateDataForRenderer( mitk::BaseRenderer *renderer )
882 LocalStorage *localStorage = m_LSH.GetLocalStorage(renderer);
884 OdfDisplayGeometry dispGeo = MeasureDisplayedGeometry( renderer);
887 if ( (localStorage->m_LastUpdateTime >= m_DataNode->GetMTime()) //was the node modified?
888 && (localStorage->m_LastUpdateTime >= m_DataNode->GetPropertyList()->GetMTime()) //was a property modified?
889 && (localStorage->m_LastUpdateTime >= m_DataNode->GetPropertyList(renderer)->GetMTime())
890 && dispGeo.Equals(m_LastDisplayGeometry))
895 localStorage->m_LastUpdateTime.Modified();
897 if(!IsVisibleOdfs(renderer))
899 localStorage->m_OdfsActors[0]->VisibilityOff();
900 localStorage->m_OdfsActors[1]->VisibilityOff();
901 localStorage->m_OdfsActors[2]->VisibilityOff();
905 localStorage->m_OdfsActors[0]->VisibilityOn();
906 localStorage->m_OdfsActors[1]->VisibilityOn();
907 localStorage->m_OdfsActors[2]->VisibilityOn();
909 m_OdfSource->SetAdditionalScale(GetMinImageSpacing(GetIndex(renderer)));
910 ApplyPropertySettings();
912 for(unsigned iter=0; iter<3; ++iter)
914 if( m_toggleColourisationMode == true )
916 localStorage-> m_OdfsMappers[iter]-> SetColorModeToDirectScalars();
917 localStorage-> m_OdfsMappers[iter]-> SelectColorArray("GLYPH_COLORS");
918 localStorage-> m_OdfsMappers[iter]-> ScalarVisibilityOn();
919 localStorage-> m_OdfsMappers[iter]-> SetScalarMode(3); // 0=Default, 1=UsePointData, 2=UseCellData, _3=UsePointFieldData_
920 // 4=UseCellFieldData, 5=UseFieldData, rot oder weiss, manchmal schwarz danach ?... double [0;1] statt char [0;255] ?
922 else if ( m_toggleColourisationMode == false )
924 localStorage-> m_OdfsMappers[iter]-> SetColorModeToDefault();
925 localStorage-> m_OdfsMappers[iter]-> SelectColorArray("vector");
926 localStorage-> m_OdfsMappers[iter]-> ScalarVisibilityOn();
927 localStorage-> m_OdfsMappers[iter]-> SetScalarMode(0); // 0=Default
931 Slice(renderer, dispGeo);
932 m_LastDisplayGeometry = dispGeo;
936 template<class T, int N>
937 double mitk::OdfVtkMapper2D<T,N>::GetMinImageSpacing( int index )
939 // Spacing adapted scaling
941 m_VtkImage->GetSpacing(spacing);
942 double min = spacing[0];
946 min = min > spacing[1] ? spacing[1] : min;
951 min = min > spacing[2] ? spacing[2] : min;
956 min = min > spacing[2] ? spacing[2] : min;
961 template<class T, int N>
962 void mitk::OdfVtkMapper2D<T,N>
963 ::ApplyPropertySettings()
965 this->GetDataNode()->GetFloatProperty( "Scaling", m_Scaling );
966 this->GetDataNode()->GetIntProperty( "ShowMaxNumber", m_ShowMaxNumber );
968 OdfNormalizationMethodProperty* nmp = dynamic_cast<OdfNormalizationMethodProperty*>(this->GetDataNode()->GetProperty( "Normalization" ));
970 m_Normalization = nmp->GetNormalization();
972 OdfScaleByProperty* sbp = dynamic_cast<OdfScaleByProperty*>(this->GetDataNode()->GetProperty( "ScaleBy" ));
974 m_ScaleBy = sbp->GetScaleBy();
976 this->GetDataNode()->GetFloatProperty( "IndexParam1", m_IndexParam1);
977 this->GetDataNode()->GetFloatProperty( "IndexParam2", m_IndexParam2);
979 this-> GetDataNode()-> GetBoolProperty( "DiffusionCore.Rendering.OdfVtkMapper.SwitchTensorView", m_toggleTensorEllipsoidView );
980 this-> GetDataNode()-> GetBoolProperty( "DiffusionCore.Rendering.OdfVtkMapper.ColourisationModeBit", m_toggleColourisationMode );
981 this-> GetDataNode()-> GetBoolProperty( "DiffusionCore.Rendering.OdfVtkMapper.RandomModeBit", m_toggleGlyphPlacementMode);
984 template <class T, int N>
985 bool mitk::OdfVtkMapper2D<T,N>
986 ::IsPlaneRotated(mitk::BaseRenderer* renderer)
988 PlaneGeometry::ConstPointer worldPlaneGeometry = renderer->GetCurrentWorldPlaneGeometry();
991 Vector3D normal = worldPlaneGeometry->GetNormal(); normal.Normalize();
992 vnl2vtk( normal.GetVnlVector(), vnormal );
994 mitk::Image* currentImage = dynamic_cast<mitk::Image* >( this->GetDataNode()->GetData() );
995 if( currentImage == NULL )
997 mitk::Vector3D imageNormal0 = currentImage->GetSlicedGeometry()->GetAxisVector(0);
998 mitk::Vector3D imageNormal1 = currentImage->GetSlicedGeometry()->GetAxisVector(1);
999 mitk::Vector3D imageNormal2 = currentImage->GetSlicedGeometry()->GetAxisVector(2);
1000 imageNormal0.Normalize();
1001 imageNormal1.Normalize();
1002 imageNormal2.Normalize();
1004 double eps = 0.000001; // Did you mean: std::numeric_limits<T>::epsilon(); ?
1006 if( fabs(fabs(dot_product(normal.GetVnlVector(),imageNormal0.GetVnlVector()))-1) > eps )
1008 if( fabs(fabs(dot_product(normal.GetVnlVector(),imageNormal1.GetVnlVector()))-1) > eps )
1010 if( fabs(fabs(dot_product(normal.GetVnlVector(),imageNormal2.GetVnlVector()))-1) > eps )
1017 template<class T, int N>
1018 void mitk::OdfVtkMapper2D<T,N>
1019 ::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* /*renderer*/, bool /*overwrite*/)
1021 node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 150 ) );
1022 node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) );
1023 node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New());
1024 node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New());
1025 node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2));
1026 node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1));
1027 node->SetProperty( "visible", mitk::BoolProperty::New( true ) );
1028 node->SetProperty( "VisibleOdfs_T", mitk::BoolProperty::New( false ) );
1029 node->SetProperty( "VisibleOdfs_C", mitk::BoolProperty::New( false ) );
1030 node->SetProperty( "VisibleOdfs_S", mitk::BoolProperty::New( false ) );
1031 node->SetProperty( "layer", mitk::IntProperty::New(100));
1032 node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) );
1033 node-> AddProperty( "DiffusionCore.Rendering.OdfVtkMapper.SwitchTensorView", mitk::BoolProperty::New( false) );
1034 node-> AddProperty( "DiffusionCore.Rendering.OdfVtkMapper.ColourisationModeBit", mitk::BoolProperty::New( false ) );
1035 node-> AddProperty( "DiffusionCore.Rendering.OdfVtkMapper.RandomModeBit", mitk::BoolProperty::New( true ) );
1038 #endif // __mitkOdfVtkMapper2D_txx__