mitk: mitkVectorImageMapper2D.cpp Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 Program:   Medical Imaging & Interaction Toolkit
00004 Language:  C++
00005 Date:      $Date$
00006 Version:   $Revision$
00007 
00008 Copyright (c) German Cancer Research Center, Division of Medical and
00009 Biological Informatics. All rights reserved.
00010 See MITKCopyright.txt or https://www.mitk.org/copyright.html for details.
00011 
00012 This software is distributed WITHOUT ANY WARRANTY; without even
00013 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00014 PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #include "mitkVectorImageMapper2D.h"
00019 
00020 //vtk related includes
00021 #include <vtkPlane.h>
00022 #include <vtkLookupTable.h>
00023 #include <vtkScalarsToColors.h>
00024 #include <vtkImageReslice.h>
00025 #include <vtkPolyData.h>
00026 #include <vtkGlyph2D.h>
00027 #include <vtkGlyphSource2D.h>
00028 #include <vtkPolyData.h>
00029 #include <vtkPolyDataSource.h>
00030 #include <vtkPolyDataMapper.h>
00031 #include <vtkPoints.h>
00032 #include <vtkCellArray.h>
00033 #include <vtkFloatArray.h>
00034 #include <vtkPointData.h>
00035 #include <vtkCellData.h>
00036 #include <vtkDataArray.h>
00037 #include <vtkMath.h>
00038 #include <vtkLinearTransform.h>
00039 #include <vtkMatrixToLinearTransform.h>
00040 #include <vtkLookupTable.h>
00041 #include <vtkScalarsToColors.h>
00042 #include <vtkTransform.h>
00043 #include <vtkImageData.h>
00044 #include <vtkDataSetWriter.h>
00045 #include <vtkMaskedGlyph3D.h>
00046 #include <vtkMaskedGlyph2D.h>
00047 #include <vtkMatrix4x4.h>
00048 #include <vtkCutter.h>
00049 #include <vtkPlane.h>
00050 #include <vtkIndent.h>
00051 #include <vtkDataObject.h>
00052 
00053 #include <fstream>
00054 
00055 //mitk related includes
00056 #include "mitkGL.h"
00057 #include "mitkBaseRenderer.h"
00058 #include "mitkColorProperty.h"
00059 #include "mitkProperties.h"
00060 #include "mitkAbstractTransformGeometry.h"
00061 #include <mitkLookupTableProperty.h>
00062 
00063 const mitk::Image * mitk::VectorImageMapper2D::GetInput( void )
00064 {
00065   if ( m_Image.IsNotNull() )
00066     return m_Image;
00067   else
00068     return dynamic_cast<const mitk::Image*>( this->GetData() );
00069 }
00070 
00071 
00072 void mitk::VectorImageMapper2D::Paint( mitk::BaseRenderer * renderer )
00073 {
00074   //std::cout << "2d vector mapping..." << std::endl;
00075   if ( IsVisible( renderer ) == false )
00076     return ;
00077 
00078   mitk::Image::Pointer input = const_cast<mitk::Image*>( this->GetInput() );
00079 
00080   if ( input.IsNull() )
00081     return ;
00082 
00083 
00084   mitk::PlaneGeometry::Pointer worldPlaneGeometry2D = dynamic_cast< mitk::PlaneGeometry*>( const_cast<mitk::Geometry2D*>( renderer->GetCurrentWorldGeometry2D() ) );
00085   assert( worldPlaneGeometry2D.IsNotNull() );
00086 
00087   vtkImageData* vtkImage = input->GetVtkImageData( this->GetCurrentTimeStep( input, renderer ) );
00088 
00089   //
00090   // set up the cutter orientation according to the current geometry of
00091   // the renderers plane
00092   //
00093   Point3D point;
00094   Vector3D normal;
00095   Geometry2D::ConstPointer worldGeometry = renderer->GetCurrentWorldGeometry2D();
00096   PlaneGeometry::ConstPointer worldPlaneGeometry = dynamic_cast<const PlaneGeometry*>( worldGeometry.GetPointer() );
00097 
00098   if ( worldPlaneGeometry.IsNotNull() )
00099   {
00100     // set up vtkPlane according to worldGeometry
00101     point = worldPlaneGeometry->GetOrigin();
00102     normal = worldPlaneGeometry->GetNormal(); normal.Normalize();
00103     m_Plane->SetTransform( (vtkAbstractTransform*)NULL );
00104   }
00105   else
00106   {
00107     itkWarningMacro( << "worldPlaneGeometry is NULL!" );
00108     return ;
00109   }
00110 
00111   vtkFloatingPointType vp[ 3 ], vp_slice[ 3 ], vnormal[ 3 ];
00112   vnl2vtk( point.Get_vnl_vector(), vp );
00113   vnl2vtk( normal.Get_vnl_vector(), vnormal );
00114   //std::cout << "Origin: " << vp[0] <<" "<< vp[1] <<" "<< vp[2] << std::endl;
00115   //std::cout << "Normal: " << vnormal[0] <<" "<< vnormal[1] <<" "<< vnormal[2] << std::endl;
00116 
00117   //normally, we would need to transform the surface and cut the transformed surface with the cutter.
00118   //This might be quite slow. Thus, the idea is, to perform an inverse transform of the plane instead.
00119   //@todo It probably does not work for scaling operations yet:scaling operations have to be
00120   //dealed with after the cut is performed by scaling the contour.
00121   vtkLinearTransform * vtktransform = GetDataNode() ->GetVtkTransform();
00122 
00123   vtkTransform* world2vtk = vtkTransform::New();
00124   world2vtk->Identity();
00125   world2vtk->Concatenate(vtktransform->GetLinearInverse());
00126   double myscale[3];
00127   world2vtk->GetScale(myscale);
00128   world2vtk->PostMultiply();
00129   world2vtk->Scale(1/myscale[0],1/myscale[1],1/myscale[2]);
00130   world2vtk->TransformPoint( vp, vp );
00131   world2vtk->TransformNormalAtPoint( vp, vnormal, vnormal );
00132   world2vtk->Delete();
00133 
00134     // vtk works in axis align coords
00135   // thus the normal also must be axis align, since 
00136   // we do not allow arbitrary cutting through volume
00137   //
00138   // vnormal should already be axis align, but in order
00139   // to get rid of precision effects, we set the two smaller
00140   // components to zero here
00141   int dims[3];
00142   vtkImage->GetDimensions(dims);
00143   double spac[3];
00144   vtkImage->GetSpacing(spac);
00145   vp_slice[0] = vp[0];
00146   vp_slice[1] = vp[1];
00147   vp_slice[2] = vp[2];
00148   if(fabs(vnormal[0]) > fabs(vnormal[1]) && fabs(vnormal[0]) > fabs(vnormal[2]) )
00149   {
00150     if(fabs(vp_slice[0]/spac[0]) < 0.4)
00151       vp_slice[0] = 0.4*spac[0];
00152     if(fabs(vp_slice[0]/spac[0]) > (dims[0]-1)-0.4)
00153       vp_slice[0] = ((dims[0]-1)-0.4)*spac[0];
00154     vnormal[1] = 0;
00155     vnormal[2] = 0;
00156   }
00157 
00158   if(fabs(vnormal[1]) > fabs(vnormal[0]) && fabs(vnormal[1]) > fabs(vnormal[2]) )
00159   {
00160     if(fabs(vp_slice[1]/spac[1]) < 0.4)
00161       vp_slice[1] = 0.4*spac[1];
00162     if(fabs(vp_slice[1]/spac[1]) > (dims[1]-1)-0.4)
00163       vp_slice[1] = ((dims[1]-1)-0.4)*spac[1];
00164     vnormal[0] = 0;
00165     vnormal[2] = 0;
00166   }
00167 
00168   if(fabs(vnormal[2]) > fabs(vnormal[1]) && fabs(vnormal[2]) > fabs(vnormal[0]) )
00169   {
00170     if(fabs(vp_slice[2]/spac[2]) < 0.4)
00171       vp_slice[2] = 0.4*spac[2];
00172     if(fabs(vp_slice[2]/spac[2]) > (dims[2]-1)-0.4)
00173       vp_slice[2] = ((dims[2]-1)-0.4)*spac[2];
00174     vnormal[0] = 0;
00175     vnormal[1] = 0;
00176   }
00177 
00178 
00179   m_Plane->SetOrigin( vp_slice );
00180   m_Plane->SetNormal( vnormal );
00181 
00182   vtkPolyData* cuttedPlane;
00183   if(!( (dims[0] == 1 && vnormal[0] != 0) || 
00184         (dims[1] == 1 && vnormal[1] != 0) ||
00185         (dims[2] == 1 && vnormal[2] != 0) ))
00186   {
00187     m_Cutter->SetCutFunction( m_Plane );
00188     m_Cutter->SetInput( vtkImage );
00189     m_Cutter->GenerateCutScalarsOff();
00190     m_Cutter->Update();
00191     cuttedPlane = m_Cutter->GetOutput();
00192   }
00193   else
00194   {
00195     // cutting of a 2D-Volume does not work, 
00196     // so we have to build up our own polydata object
00197     cuttedPlane = vtkPolyData::New();
00198     vtkPoints* points = vtkPoints::New();
00199     points->SetNumberOfPoints(vtkImage->GetNumberOfPoints());
00200     for(int i=0; i<vtkImage->GetNumberOfPoints(); i++)
00201       points->SetPoint(i, vtkImage->GetPoint(i));
00202     cuttedPlane->SetPoints(points);
00203     vtkFloatArray* pointdata = vtkFloatArray::New();
00204     int comps  = vtkImage->GetPointData()->GetScalars()->GetNumberOfComponents();
00205     pointdata->SetNumberOfComponents(comps);
00206     int tuples = vtkImage->GetPointData()->GetScalars()->GetNumberOfTuples();
00207     pointdata->SetNumberOfTuples(tuples);
00208     for(int i=0; i<tuples; i++)
00209       pointdata->SetTuple(i,vtkImage->GetPointData()->GetScalars()->GetTuple(i));
00210     pointdata->SetName( "vector" );
00211     cuttedPlane->GetPointData()->AddArray(pointdata);
00212   }
00213 
00214   if ( cuttedPlane->GetNumberOfPoints() != 0)
00215   {
00216     //
00217     // make sure, that we have point data with more than 1 component (as vectors)
00218     //
00219     vtkPointData * pointData = cuttedPlane->GetPointData();
00220     if ( pointData == NULL )
00221     {
00222       itkWarningMacro( << "no point data associated with cutters result!" );
00223       return ;
00224     }
00225     if ( pointData->GetNumberOfArrays() == 0 )
00226     {
00227       itkWarningMacro( << "point data returned by cutter doesn't have any arrays associated!" );
00228       return ;
00229     }
00230     else if ( pointData->GetArray(0)->GetNumberOfComponents() <= 1)
00231     {
00232       itkWarningMacro( << "number of components <= 1!" );
00233       return;
00234     }
00235     else if ( pointData->GetArrayName( 0 ) == NULL )
00236     {
00237       pointData->GetArray( 0 ) ->SetName( "vector" );
00238       //std::cout << "array name = vectors now" << std::endl;
00239     }
00240     //std::cout << "  projecting..."<< std::endl;
00241 
00242     //
00243     // constrain the vectors to lie on the plane, which means to remove the vector component,
00244     // which is orthogonal to the plane.
00245     //
00246     vtkIdType numPoints, pointId;
00247     numPoints = cuttedPlane->GetNumberOfPoints();
00248     vtkDataArray* inVectors = cuttedPlane->GetPointData()->GetVectors( "vector" );
00249     assert( inVectors != NULL );
00250     vtkFloatArray* vectorMagnitudes = vtkFloatArray::New();
00251     vectorMagnitudes->SetName("vectorMagnitudes");
00252     vectorMagnitudes->SetNumberOfComponents(1);
00253     vectorMagnitudes->SetNumberOfValues(numPoints);
00254     vectorMagnitudes->SetNumberOfTuples(numPoints);
00255     vtkFloatingPointType inVector[ 3 ], outVector[3], wnormal[3]; //, tmpVector[ 3 ], outVector[ 3 ];
00256     vtkFloatingPointType k = 0.0;
00257     vnl2vtk( normal.Get_vnl_vector(), wnormal );
00258     vtkMath::Normalize( wnormal );
00259     bool normalizeVecs;
00260     m_DataNode->GetBoolProperty( "NormalizeVecs", normalizeVecs );
00261     for ( pointId = 0; pointId < numPoints; ++pointId )
00262     {
00263       inVectors->GetTuple( pointId, inVector );
00264       if(normalizeVecs)
00265       {
00266         vnl_vector<double> tmp(3);
00267         vtk2vnl(inVector, tmp);
00268         tmp.normalize();
00269         vnl2vtk(tmp, inVector);
00270       }
00271       k = vtkMath::Dot( wnormal, inVector );
00272       // Remove non orthogonal component.
00273       outVector[ 0 ] = inVector[ 0 ] - ( wnormal[ 0 ] * k );
00274       outVector[ 1 ] = inVector[ 1 ] - ( wnormal[ 1 ] * k );
00275       outVector[ 2 ] = inVector[ 2 ] - ( wnormal[ 2 ] * k );
00276       inVectors->SetTuple( pointId, outVector );
00277 
00278       // ?? this was set to norm(inVector) before, but outVector made more sense to me
00279       vectorMagnitudes->SetValue( pointId, vtkMath::Norm( outVector ) );
00280 
00281 
00282       //std::cout << "method old: " << inVector[0] <<", " << inVector[1] << ", "<<inVector[2] << ", method new: " << outVector[0] << ", "<< outVector[1] << ", "<< outVector[2] << std::endl;
00283     }
00284     pointData->AddArray(vectorMagnitudes);
00285     pointData->CopyAllOn();
00286 
00287     //pointData->PrintSelf(std::cout, vtkIndent(4));
00288     //std::cout << "  ...done!"<< std::endl;
00289     //std::cout << "  glyphing..."<< std::endl;
00290 
00291     // call glyph2D to generate 2D glyphs for each of the
00292     // vectors
00293     vtkGlyphSource2D* glyphSource = vtkGlyphSource2D::New();
00294     //glyphSource->SetGlyphTypeToDash();
00295     glyphSource->DashOn();
00296     //glyphSource->SetScale( 0.1 );
00297     //glyphSource->SetScale2( .5 );
00298     //glyphSource->SetCenter( 0.5, 0.5, 0.5 );
00299     glyphSource->CrossOff();
00300     //glyphSource->FilledOff();
00301     //glyphSource->Update();
00302 
00303     double spacing[3];
00304     vtkImage->GetSpacing(spacing);
00305     double min = spacing[0];
00306     min = min > spacing[1] ? spacing[1] : min;
00307     min = min > spacing[2] ? spacing[2] : min;
00308 
00309     float scale = 1;
00310     mitk::FloatProperty::Pointer mitkScaleProp = dynamic_cast<mitk::FloatProperty*>(GetDataNode()->GetProperty("Scale"));
00311     if (mitkScaleProp.IsNotNull())
00312     {
00313       scale = mitkScaleProp->GetValue();
00314     }
00315 
00316     vtkMaskedGlyph3D* glyphGenerator = vtkMaskedGlyph3D::New();
00317     glyphGenerator->SetSource( glyphSource->GetOutput() );
00318     glyphGenerator->SetInputConnection(cuttedPlane->GetProducerPort());
00319     glyphGenerator->SetInputArrayToProcess (1, 0,0, vtkDataObject::FIELD_ASSOCIATION_POINTS , "vector");
00320     glyphGenerator->SetVectorModeToUseVector();
00321     glyphGenerator->OrientOn();
00322     glyphGenerator->SetScaleFactor( min*scale );
00323     glyphGenerator->SetUseMaskPoints( true );
00324     glyphGenerator->SetRandomMode( true );
00325     glyphGenerator->SetMaximumNumberOfPoints( 128*128 );
00326 
00327     glyphGenerator->Update();
00328 
00329     vtkLookupTable* vtkLut = NULL;
00330     mitk::LookupTableProperty::Pointer mitkLutProp = dynamic_cast<mitk::LookupTableProperty*>(GetDataNode()->GetProperty("LookupTable"));
00331     if (mitkLutProp.IsNotNull())
00332     {
00333       vtkLut = mitkLutProp->GetLookupTable()->GetVtkLookupTable();
00334     }
00335 
00336     mitk::Color color;
00337     mitk::ColorProperty::Pointer mitkColorProp = dynamic_cast<mitk::ColorProperty*>(GetDataNode()->GetProperty("color"));
00338     if (mitkColorProp.IsNotNull())
00339     {
00340       color = mitkColorProp->GetColor();
00341     }
00342     else
00343     {
00344       color.SetRed(0);
00345       color.SetBlue(1);
00346       color.SetGreen(0);
00347     }
00348 
00349     float lwidth = 1;
00350     mitk::FloatProperty::Pointer mitkLWidthProp = dynamic_cast<mitk::FloatProperty*>(GetDataNode()->GetProperty("LineWidth"));
00351     if (mitkLWidthProp.IsNotNull())
00352     {
00353       lwidth = mitkLWidthProp->GetValue();
00354     }
00355 
00356     vtkTransform* trafo = vtkTransform::New();
00357     trafo->Identity();
00358     trafo->Concatenate(vtktransform);
00359     trafo->PreMultiply();
00360     double myscale[3];
00361     trafo->GetScale(myscale);
00362     trafo->Scale(1/myscale[0],1/myscale[1],1/myscale[2]);
00363     
00364     this->PaintCells( glyphGenerator->GetOutput(), renderer->GetCurrentWorldGeometry2D(), renderer->GetDisplayGeometry(), trafo, renderer, NULL/*vtkLut*/, color, lwidth, spacing );
00365 
00366     vectorMagnitudes->Delete();
00367     glyphSource->Delete();
00368     glyphGenerator->Delete();
00369     trafo->Delete();
00370   }
00371   else
00372   {
00373     std::cout << "  no points cutted!"<< std::endl;
00374   }
00375   //std::cout << "...done!" << std::endl;
00376 }
00377 
00378 
00379 
00380 
00381 
00382 
00383 
00384 void mitk::VectorImageMapper2D::PaintCells( vtkPolyData* glyphs, const Geometry2D* worldGeometry, const DisplayGeometry* displayGeometry, vtkLinearTransform* vtktransform, mitk::BaseRenderer*  /*renderer*/, vtkScalarsToColors *lut, mitk::Color color, float lwidth, vtkFloatingPointType *spacing )
00385 {
00386 
00387   vtkPoints * points = glyphs->GetPoints();
00388   vtkPointData * vpointdata = glyphs->GetPointData();
00389   vtkDataArray* vpointscalars = vpointdata->GetArray("vectorMagnitudes");
00390   //vtkDataArray* vpointpositions = vpointdata->GetArray("pointPositions");
00391   assert(vpointscalars != NULL);
00392   //std::cout << "  Scalars range 2d:" << vpointscalars->GetRange()[0] << " " << vpointscalars->GetRange()[0] << std::endl;
00393 
00394   Point3D p;
00395   Point2D p2d;
00396   vtkIdList* idList;
00397   vtkCell* cell;
00398 
00399   vtkFloatingPointType offset[3]; 
00400   for (unsigned int i = 0; i < 3; ++i) 
00401   {
00402     offset[i] = 0; 
00403   }
00404 
00405   vtkIdType numCells = glyphs->GetNumberOfCells();
00406   for ( vtkIdType cellId = 0; cellId < numCells; ++cellId )
00407   {
00408     vtkFloatingPointType vp[ 3 ];
00409 
00410     cell = glyphs->GetCell( cellId );
00411     idList = cell->GetPointIds();
00412 
00413     int numPoints = idList->GetNumberOfIds();
00414 
00415     if(numPoints == 1)
00416     {
00417       //take transformation via vtktransform into account
00418       vtkFloatingPointType pos[ 3 ],vp_raster[3];
00419       points->GetPoint( idList->GetId( 0 ), vp );
00420       vp_raster[0] = vtkMath::Round(vp[0]/spacing[0])*spacing[0];
00421       vp_raster[1] = vtkMath::Round(vp[1]/spacing[1])*spacing[1];
00422       vp_raster[2] = vtkMath::Round(vp[2]/spacing[2])*spacing[2];
00423       vtktransform->TransformPoint( vp_raster, pos );
00424       offset[0] = pos[0] - vp[0];
00425       offset[1] = pos[1] - vp[1];
00426       offset[2] = pos[2] - vp[2];
00427     }
00428     else
00429     {
00430       glLineWidth(lwidth);
00431       glBegin ( GL_LINE_LOOP );
00432       
00433       for ( int pointNr = 0; pointNr < numPoints ;++pointNr )
00434       {
00435         points->GetPoint( idList->GetId( pointNr ), vp );
00436 
00437         vp[0] = vp[0] + offset[0];
00438         vp[1] = vp[1] + offset[1];
00439         vp[2] = vp[2] + offset[2];
00440 
00441         vtkFloatingPointType tmp[ 3 ];
00442         vtktransform->TransformPoint( vp,tmp );
00443 
00444         vtk2itk( vp, p );
00445 
00446         //convert 3D point (in mm) to 2D point on slice (also in mm)
00447         worldGeometry->Map( p, p2d );
00448 
00449         //convert point (until now mm and in worldcoordinates) to display coordinates (units )
00450         displayGeometry->WorldToDisplay( p2d, p2d );
00451 
00452         if ( lut != NULL )
00453         {
00454           // color each point according to point data
00455           vtkFloatingPointType * color;
00456 
00457           if ( vpointscalars != NULL )
00458           {
00459             vpointscalars->GetComponent( pointNr, 0 );
00460             color = lut->GetColor( vpointscalars->GetComponent( idList->GetId( pointNr ), 0 ) );
00461             glColor3f( color[ 0 ], color[ 1 ], color[ 2 ] );
00462           }
00463         }
00464         else
00465         {
00466           glColor3f( color.GetRed(), color.GetGreen(), color.GetBlue() );
00467         }
00468 
00469         //std::cout <<  idList->GetId( pointNr )<< ": " << p2d[0]<< " "<< p2d[1] << std::endl;
00470         //draw the line
00471         glVertex2f( p2d[ 0 ], p2d[ 1 ] );
00472 
00473 
00474 
00475       }
00476       glEnd ();
00477     }
00478   }
00479 }
00480 
00481 
00482 mitk::VectorImageMapper2D::VectorImageMapper2D()
00483 {
00484   m_LUT = NULL;
00485   m_Plane = vtkPlane::New();
00486   m_Cutter = vtkCutter::New();
00487 
00488   m_Cutter->SetCutFunction( m_Plane );
00489   m_Cutter->GenerateValues( 1, 0, 1 );
00490 }
00491 
00492 
00493 mitk::VectorImageMapper2D::~VectorImageMapper2D()
00494 {
00495   if ( m_LUT != NULL )
00496     m_LUT->Delete();
00497   if ( m_Plane != NULL )
00498     m_Plane->Delete();
00499   if ( m_Cutter != NULL )
00500     m_Cutter->Delete();
00501 }
00502 
00503 int mitk::VectorImageMapper2D::GetCurrentTimeStep( mitk::BaseData* data, mitk::BaseRenderer* renderer )
00504 {
00505   //
00506   // get the TimeSlicedGeometry of the input object
00507   //
00508   const TimeSlicedGeometry * dataTimeGeometry = data->GetUpdatedTimeSlicedGeometry();
00509   if ( ( dataTimeGeometry == NULL ) || ( dataTimeGeometry->GetTimeSteps() == 0 ) )
00510   {
00511     itkWarningMacro( << "geometry of the given data object isn't a mitk::TimeSlicedGeometry, or the number of time steps is 0!" );
00512     return 0;
00513   }
00514 
00515   //
00516   // get the world time
00517   //
00518   Geometry2D::ConstPointer worldGeometry = renderer->GetCurrentWorldGeometry2D();
00519   assert( worldGeometry.IsNotNull() );
00520   ScalarType time = worldGeometry->GetTimeBounds() [ 0 ];
00521 
00522   //
00523   // convert the world time to time steps of the input object
00524   //
00525   int timestep = 0;
00526   if ( time > ScalarTypeNumericTraits::NonpositiveMin() )
00527     timestep = dataTimeGeometry->MSToTimeStep( time );
00528   if ( dataTimeGeometry->IsValidTime( timestep ) == false )
00529   {
00530     itkWarningMacro( << timestep << " is not a valid time of the given data object!" );
00531     return 0;
00532   }
00533   return timestep;
00534 }
00535