mitk: mitkExtrudedContour.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 
00019 #include "mitkExtrudedContour.h"
00020 #include "mitkVector.h"
00021 #include "mitkBaseProcess.h"
00022 
00023 #include <itkSmartPointerForwardReference.txx>
00024 
00025 #include <vtkLinearExtrusionFilter.h>
00026 #include <vtkPolyData.h>
00027 #include <vtkPoints.h>
00028 #include <vtkCellArray.h>
00029 #include <vtkPlanes.h>
00030 #include <vtkClipPolyData.h>
00031 
00032 #include <vtkPolygon.h>
00033 #include <vtkFloatArray.h>
00034 #include <vtkDoubleArray.h>
00035 #include <vtkPlane.h>
00036 
00037 
00038 //vtkButterflySubdivisionFilter * subdivs;
00039 #include <vtkSampleFunction.h>
00040 #include <vtkTriangleFilter.h>
00041 #include <vtkLinearSubdivisionFilter.h>
00042 #include <vtkImageData.h>
00043 #include <vtkCubeSource.h>
00044 
00045 mitk::ExtrudedContour::ExtrudedContour() 
00046   : m_Contour(NULL), m_ClippingGeometry(NULL), m_AutomaticVectorGeneration(false)
00047 {
00048   GetTimeSlicedGeometry()->Initialize(1);
00049 
00050   FillVector3D(m_Vector, 0.0, 0.0, 1.0);
00051   m_RightVector.Fill(0.0);
00052 
00053   m_ExtrusionFilter = vtkLinearExtrusionFilter::New();
00054 
00055   m_ExtrusionFilter->CappingOff();
00056   m_ExtrusionFilter->SetExtrusionTypeToVectorExtrusion();
00057 
00058 #if ((VTK_MAJOR_VERSION > 4) || ((VTK_MAJOR_VERSION==4) && (VTK_MINOR_VERSION>=4) ))
00059   double vtkvector[3]={0,0,1};
00060 #else
00061   float vtkvector[3]={0,0,1};
00062 #endif
00063   // set extrusion vector
00064   m_ExtrusionFilter->SetVector(vtkvector);
00065 
00066   m_TriangleFilter = vtkTriangleFilter::New();
00067   m_TriangleFilter->SetInput(m_ExtrusionFilter->GetOutput());
00068 
00069   m_SubdivisionFilter = vtkLinearSubdivisionFilter::New();
00070   m_SubdivisionFilter->SetInput(m_TriangleFilter->GetOutput());
00071   m_SubdivisionFilter->SetNumberOfSubdivisions(4);
00072 
00073   m_ClippingBox = vtkPlanes::New();
00074   m_ClipPolyDataFilter = vtkClipPolyData::New();
00075   m_ClipPolyDataFilter->SetInput(m_SubdivisionFilter->GetOutput());
00076   m_ClipPolyDataFilter->SetClipFunction(m_ClippingBox);
00077   m_ClipPolyDataFilter->InsideOutOn();
00078 
00079   m_Polygon = vtkPolygon::New();
00080   
00081   m_ProjectionPlane = mitk::PlaneGeometry::New();
00082 }
00083 
00084 mitk::ExtrudedContour::~ExtrudedContour()
00085 {
00086   m_ClipPolyDataFilter->Delete();
00087   m_ClippingBox->Delete();
00088   m_SubdivisionFilter->Delete();
00089   m_TriangleFilter->Delete();
00090   m_ExtrusionFilter->Delete();
00091   m_Polygon->Delete();
00092 }
00093 
00094 bool mitk::ExtrudedContour::IsInside(const Point3D& worldPoint) const
00095 {
00096 #if ((VTK_MAJOR_VERSION > 4) || ((VTK_MAJOR_VERSION==4) && (VTK_MINOR_VERSION>=4) ))
00097   static double polygonNormal[3]={0.0,0.0,1.0};
00098 #else
00099   static float polygonNormal[3]={0.0,0.0,1.0};
00100 #endif
00101 
00102 
00103   // project point onto plane
00104   float xt[3];
00105   itk2vtk(worldPoint, xt);
00106 
00107   xt[0] = worldPoint[0]-m_Origin[0];
00108   xt[1] = worldPoint[1]-m_Origin[1];
00109   xt[2] = worldPoint[2]-m_Origin[2];
00110 
00111   float dist=xt[0]*m_Normal[0]+xt[1]*m_Normal[1]+xt[2]*m_Normal[2];
00112   xt[0] -= dist*m_Normal[0];
00113   xt[1] -= dist*m_Normal[1];
00114   xt[2] -= dist*m_Normal[2];
00115 
00116 #if ((VTK_MAJOR_VERSION > 4) || ((VTK_MAJOR_VERSION==4) && (VTK_MINOR_VERSION>=4) ))
00117     double x[3];
00118 #else
00119    float x[3];
00120 #endif
00121 
00122   x[0] = xt[0]*m_Right[0]+xt[1]*m_Right[1]+xt[2]*m_Right[2];
00123   x[1] = xt[0]*m_Down[0] +xt[1]*m_Down[1] +xt[2]*m_Down[2];
00124   x[2] = 0;
00125 
00126 #if ((VTK_MAJOR_VERSION > 4) || ((VTK_MAJOR_VERSION==4) && (VTK_MINOR_VERSION>=4) ))
00127   // determine whether it's in the selection loop and then evaluate point
00128   // in polygon only if absolutely necessary.
00129   if ( x[0] >= this->m_ProjectedContourBounds[0] && x[0] <= this->m_ProjectedContourBounds[1] &&
00130     x[1] >= this->m_ProjectedContourBounds[2] && x[1] <= this->m_ProjectedContourBounds[3] &&
00131     this->m_Polygon->PointInPolygon(x, m_Polygon->Points->GetNumberOfPoints(),
00132     ((vtkDoubleArray *)this->m_Polygon->Points->GetData())->GetPointer(0), 
00133     (double*)const_cast<mitk::ExtrudedContour*>(this)->m_ProjectedContourBounds, polygonNormal) == 1 )
00134     return true;
00135   else
00136     return false;
00137 #else
00138   // determine whether it's in the selection loop and then evaluate point
00139   // in polygon only if absolutely necessary.
00140   if ( x[0] >= this->m_ProjectedContourBounds[0] && x[0] <= this->m_ProjectedContourBounds[1] &&
00141     x[1] >= this->m_ProjectedContourBounds[2] && x[1] <= this->m_ProjectedContourBounds[3] &&
00142     this->m_Polygon->PointInPolygon(x, m_Polygon->Points->GetNumberOfPoints(),
00143     ((vtkFloatArray *)this->m_Polygon->Points->GetData())->GetPointer(0), 
00144     const_cast<mitk::ExtrudedContour*>(this)->m_ProjectedContourBounds, polygonNormal) == 1 )
00145     return true;
00146   else
00147     return false;
00148 #endif
00149 
00150 }
00151 
00152 mitk::ScalarType mitk::ExtrudedContour::GetVolume()
00153 {
00154   return -1.0;
00155 }
00156 
00157 void mitk::ExtrudedContour::UpdateOutputInformation()
00158 {
00159   if ( this->GetSource() )
00160   {
00161     this->GetSource()->UpdateOutputInformation();
00162   }
00163   if(GetMTime() > m_LastCalculateExtrusionTime)
00164   {
00165     BuildGeometry();
00166     BuildSurface();
00167   }
00168   //if ( ( m_CalculateBoundingBox ) && ( m_PolyDataSeries.size() > 0 ) )
00169   //  CalculateBoundingBox();
00170 }
00171 
00172 void mitk::ExtrudedContour::BuildSurface()
00173 {
00174   if(m_Contour.IsNull())
00175   {
00176     SetVtkPolyData(NULL);
00177     return;
00178   }
00179 
00180   // set extrusion contour
00181   vtkPolyData *polyData = vtkPolyData::New();
00182   vtkCellArray *polys   = vtkCellArray::New();
00183 
00184   polys->InsertNextCell(m_Polygon->GetPointIds());
00185   polyData->SetPoints(m_Polygon->GetPoints());
00186 
00187   //float vtkpoint[3];
00188   //unsigned int i, numPts = m_Polygon->GetNumberOfPoints();
00189   //for(i=0; i<numPts; ++i)
00190   //{
00191   //  float * vtkpoint = this->m_Polygon->Points->GetPoint(i);
00192 
00193   //  pointids[i]=loopPoints->InsertNextPoint(vtkpoint);
00194   //}
00195   //polys->InsertNextCell( i, pointids );
00196   //delete [] pointids;
00197 
00198   //polyData->SetPoints( loopPoints );
00199   polyData->SetPolys( polys );
00200   polys->Delete();
00201 
00202   m_ExtrusionFilter->SetInput(polyData);
00203 
00204   polyData->Delete();
00205 
00206   // set extrusion scale factor
00207   m_ExtrusionFilter->SetScaleFactor(GetGeometry()->GetExtentInMM(2));
00208   SetVtkPolyData(m_SubdivisionFilter->GetOutput());
00209   //if(m_ClippingGeometry.IsNull())
00210   //{
00211   //  SetVtkPolyData(m_SubdivisionFilter->GetOutput());
00212   //}
00213   //else
00214   //{
00215   //  m_ClipPolyDataFilter->SetInput(m_SubdivisionFilter->GetOutput());
00216   //  mitk::BoundingBox::BoundsArrayType bounds=m_ClippingGeometry->GetBounds();
00217 
00218   //  m_ClippingBox->SetBounds(bounds[0], bounds[1], bounds[2], bounds[3], bounds[4], bounds[5]);
00219   //  m_ClippingBox->SetTransform(GetGeometry()->GetVtkTransform());
00220   //  m_ClipPolyDataFilter->SetClipFunction(m_ClippingBox);
00221   //  m_ClipPolyDataFilter->SetValue(0);
00222   //  SetVtkPolyData(m_ClipPolyDataFilter->GetOutput());
00223   //}
00224 
00225   m_LastCalculateExtrusionTime.Modified();
00226 }
00227 
00228 void mitk::ExtrudedContour::BuildGeometry()
00229 {
00230   if(m_Contour.IsNull())
00231     return;
00232 
00233 //  Initialize(1);
00234 
00235   Vector3D nullvector; nullvector.Fill(0.0);
00236   float xProj[3];
00237   unsigned int i;
00238   unsigned int numPts = 20; //m_Contour->GetNumberOfPoints();
00239   mitk::Contour::PathPointer path = m_Contour->GetContourPath();
00240   mitk::Contour::PathType::InputType cstart = path->StartOfInput();
00241   mitk::Contour::PathType::InputType cend   = path->EndOfInput();
00242   mitk::Contour::PathType::InputType cstep  = (cend-cstart)/numPts;
00243   mitk::Contour::PathType::InputType ccur;
00244 
00245   // Part I: guarantee/calculate legal vectors
00246 
00247   m_Vector.Normalize();
00248   itk2vtk(m_Vector, m_Normal);
00249   // check m_Vector
00250   if(mitk::Equal(m_Vector, nullvector) || m_AutomaticVectorGeneration)
00251   {
00252     if ( m_AutomaticVectorGeneration == false)
00253       itkWarningMacro("Extrusion vector is 0 ("<< m_Vector << "); trying to use normal of polygon");
00254 
00255     vtkPoints *loopPoints = vtkPoints::New();
00256     //mitk::Contour::PointsContainerIterator pointsIt = m_Contour->GetPoints()->Begin();
00257 #if ((VTK_MAJOR_VERSION > 4) || ((VTK_MAJOR_VERSION==4) && (VTK_MINOR_VERSION>=4) ))
00258     double vtkpoint[3];
00259 #else
00260     float vtkpoint[3];
00261 #endif
00262     unsigned int i=0;
00263     for(i=0, ccur=cstart; i<numPts; ++i, ccur+=cstep)
00264     {
00265       itk2vtk(path->Evaluate(ccur), vtkpoint);
00266       loopPoints->InsertNextPoint(vtkpoint);
00267     }
00268 
00269     // Make sure points define a loop with a m_Normal
00270     vtkPolygon::ComputeNormal(loopPoints, m_Normal);
00271     loopPoints->Delete();
00272 
00273     vtk2itk(m_Normal, m_Vector);
00274     if(mitk::Equal(m_Vector, nullvector))
00275     {
00276       itkExceptionMacro("Cannot calculate normal of polygon");
00277     }
00278   }
00279   // check m_RightVector
00280   if((mitk::Equal(m_RightVector, nullvector)) || (mitk::Equal(m_RightVector*m_Vector, 0.0)==false))
00281   {
00282     if(mitk::Equal(m_RightVector, nullvector))
00283     {
00284       itkDebugMacro("Right vector is 0. Calculating.");
00285     }
00286     else
00287     {
00288       itkWarningMacro("Right vector ("<<m_RightVector<<") not perpendicular to extrusion vector "<<m_Vector<<": "<<m_RightVector*m_Vector);
00289     }
00290     // calculate a legal m_RightVector
00291     if( mitk::Equal( m_Vector[1], 0.0f ) == false )
00292     {
00293       FillVector3D( m_RightVector, 1.0f, -m_Vector[0]/m_Vector[1], 0.0f ); 
00294       m_RightVector.Normalize();
00295     }
00296     else
00297     {
00298       FillVector3D( m_RightVector, 0.0f, 1.0f, 0.0f );
00299     }
00300   }
00301  
00302   // calculate down-vector
00303   VnlVector rightDV = m_RightVector.Get_vnl_vector(); rightDV.normalize(); vnl2vtk(rightDV, m_Right);
00304   VnlVector downDV  = vnl_cross_3d( m_Vector.Get_vnl_vector(), rightDV ); downDV.normalize();  vnl2vtk(downDV,  m_Down);
00305   
00306   // Part II: calculate plane as base for extrusion, project the contour
00307   // on this plane and store as polygon for IsInside test and BoundingBox calculation
00308 
00309   // initialize m_ProjectionPlane, yet with origin at 0
00310   m_ProjectionPlane->InitializeStandardPlane(rightDV, downDV);
00311 
00312   // create vtkPolygon from contour and simultaneously determine 2D bounds of 
00313   // contour projected on m_ProjectionPlane
00314   //mitk::Contour::PointsContainerIterator pointsIt = m_Contour->GetPoints()->Begin();
00315   m_Polygon->Points->Reset();
00316   m_Polygon->Points->SetNumberOfPoints(numPts);
00317   m_Polygon->PointIds->Reset();
00318   m_Polygon->PointIds->SetNumberOfIds(numPts);
00319   mitk::Point2D pt2d;
00320   mitk::Point3D pt3d;
00321   mitk::Point2D min, max;
00322   min.Fill(ScalarTypeNumericTraits::max());
00323   max.Fill(ScalarTypeNumericTraits::min());
00324   xProj[2]=0.0;
00325   for(i=0, ccur=cstart; i<numPts; ++i, ccur+=cstep)
00326   {
00327     pt3d.CastFrom(path->Evaluate(ccur));
00328     m_ProjectionPlane->Map(pt3d, pt2d);
00329     xProj[0]=pt2d[0];
00330     if(pt2d[0]<min[0]) min[0]=pt2d[0];
00331     if(pt2d[0]>max[0]) max[0]=pt2d[0];
00332     xProj[1]=pt2d[1];
00333     if(pt2d[1]<min[1]) min[1]=pt2d[1];
00334     if(pt2d[1]>max[1]) max[1]=pt2d[1];
00335     m_Polygon->Points->SetPoint(i, xProj);
00336     m_Polygon->PointIds->SetId(i, i);
00337   }
00338   // shift parametric origin to (0,0)
00339   for(i=0; i<numPts; ++i)
00340   {
00341 #if ((VTK_MAJOR_VERSION > 4) || ((VTK_MAJOR_VERSION==4) && (VTK_MINOR_VERSION>=4) ))
00342     double * pt = this->m_Polygon->Points->GetPoint(i);
00343 #else
00344     float * pt = this->m_Polygon->Points->GetPoint(i);
00345 #endif
00346     pt[0]-=min[0]; pt[1]-=min[1];
00347     itkDebugMacro( << i << ": (" << pt[0] << "," << pt[1] << "," << pt[2] << ")" );
00348   }
00349   this->m_Polygon->GetBounds(m_ProjectedContourBounds);
00350   //m_ProjectedContourBounds[4]=-1.0; m_ProjectedContourBounds[5]=1.0;
00351 
00352   // calculate origin (except translation along the normal) and bounds 
00353   // of m_ProjectionPlane:
00354   //  origin is composed of the minimum x-/y-coordinates of the polygon,
00355   //  bounds from the extent of the polygon, both after projecting on the plane
00356   mitk::Point3D origin;
00357   m_ProjectionPlane->Map(min, origin);
00358   ScalarType bounds[6]={0, max[0]-min[0], 0, max[1]-min[1], 0, 1};
00359   m_ProjectionPlane->SetBounds(bounds);
00360   m_ProjectionPlane->SetOrigin(origin);
00361 
00362   // Part III: initialize geometry
00363   if(m_ClippingGeometry.IsNotNull())
00364   {
00365     ScalarType min_dist=ScalarTypeNumericTraits::max(), max_dist=ScalarTypeNumericTraits::min(), dist;
00366     unsigned char i;
00367     for(i=0; i<8; ++i)
00368     {
00369       dist = m_ProjectionPlane->SignedDistance( m_ClippingGeometry->GetCornerPoint(i) );
00370       if(dist<min_dist) min_dist=dist;
00371       if(dist>max_dist) max_dist=dist;
00372     }
00373     //incorporate translation along the normal into origin
00374     origin = origin+m_Vector*min_dist;
00375     m_ProjectionPlane->SetOrigin(origin);
00376     bounds[5]=max_dist-min_dist;
00377   }
00378   else
00379     bounds[5]=20;
00380 
00381   itk2vtk(origin, m_Origin);
00382 
00383   mitk::TimeSlicedGeometry::Pointer timeGeometry = this->GetTimeSlicedGeometry();
00384   mitk::Geometry3D::Pointer g3d = timeGeometry->GetGeometry3D( 0 );
00385   assert( g3d.IsNotNull() );
00386   g3d->SetBounds(bounds);
00387   g3d->SetIndexToWorldTransform(m_ProjectionPlane->GetIndexToWorldTransform());
00388   g3d->TransferItkToVtkTransform();
00389   timeGeometry->InitializeEvenlyTimed(g3d, 1);
00390 }
00391 
00392 unsigned long mitk::ExtrudedContour::GetMTime() const
00393 {
00394   unsigned long latestTime = Superclass::GetMTime();
00395   if(m_Contour.IsNotNull())
00396   {
00397     unsigned long localTime;
00398     localTime = m_Contour->GetMTime();
00399     if(localTime > latestTime) latestTime = localTime;
00400   }
00401   return latestTime;
00402 }