mitk: mitkGeometryClipImageFilter.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 "mitkGeometryClipImageFilter.h"
00020 #include "mitkImageTimeSelector.h"
00021 #include "mitkTimeHelper.h"
00022 #include "mitkProperties.h"
00023 
00024 #include "mitkImageToItk.h"
00025 
00026 #include "itkImageRegionConstIterator.h"
00027 #include "itkImageRegionIteratorWithIndex.h"
00028 
00029 #include <limits>
00030 
00031 mitk::GeometryClipImageFilter::GeometryClipImageFilter() 
00032   : m_ClippingGeometry(NULL), 
00033     m_ClipPartAboveGeometry(true), 
00034     m_OutsideValue(0), 
00035     m_AutoOutsideValue(false), 
00036     m_LabelBothSides(false),
00037     m_AutoOrientLabels(false),
00038     m_AboveGeometryLabel(1), 
00039     m_BelowGeometryLabel(2)
00040 {
00041   this->SetNumberOfInputs(2);
00042   this->SetNumberOfRequiredInputs(2);
00043   m_InputTimeSelector  = mitk::ImageTimeSelector::New();
00044   m_OutputTimeSelector = mitk::ImageTimeSelector::New();
00045   m_ClippingGeometryData = mitk::GeometryData::New();
00046 }
00047 
00048 mitk::GeometryClipImageFilter::~GeometryClipImageFilter()
00049 {
00050 
00051 }
00052 
00053 void mitk::GeometryClipImageFilter::SetClippingGeometry(const mitk::Geometry3D* aClippingGeometry)
00054 {
00055   if(aClippingGeometry != m_ClippingGeometry.GetPointer())
00056   {
00057     m_ClippingGeometry = aClippingGeometry;
00058     m_TimeSlicedClippingGeometry = dynamic_cast<const TimeSlicedGeometry*>(aClippingGeometry);
00059     m_ClippingGeometryData->SetGeometry(const_cast<mitk::Geometry3D*>(aClippingGeometry));
00060     SetNthInput(1, m_ClippingGeometryData);
00061     Modified();
00062   }
00063 }
00064 
00065 const mitk::Geometry3D* mitk::GeometryClipImageFilter::GetClippingGeometry() const 
00066 {
00067   return m_ClippingGeometry;
00068 }
00069 
00070 void mitk::GeometryClipImageFilter::GenerateInputRequestedRegion()
00071 {
00072   Superclass::GenerateInputRequestedRegion();
00073 
00074   mitk::Image* output = this->GetOutput();
00075   mitk::Image* input = const_cast< mitk::Image * > ( this->GetInput() );
00076   if((output->IsInitialized()==false) || (m_ClippingGeometry.IsNull()))
00077     return;
00078 
00079   input->SetRequestedRegionToLargestPossibleRegion();
00080 
00081   GenerateTimeInInputRegion(output, input);
00082 }
00083 
00084 void mitk::GeometryClipImageFilter::GenerateOutputInformation()
00085 {
00086   mitk::Image::ConstPointer input = this->GetInput();
00087   mitk::Image::Pointer output = this->GetOutput();
00088 
00089   if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime()))
00090     return;
00091 
00092   itkDebugMacro(<<"GenerateOutputInformation()");
00093 
00094   unsigned int i;
00095   unsigned int *tmpDimensions = new unsigned int[input->GetDimension()];
00096 
00097   for(i=0;i<input->GetDimension();++i)
00098     tmpDimensions[i]=input->GetDimension(i);
00099 
00100   output->Initialize(input->GetPixelType(),
00101     input->GetDimension(),
00102     tmpDimensions,
00103     input->GetNumberOfChannels());
00104 
00105   delete [] tmpDimensions;
00106 
00107   output->SetGeometry(static_cast<mitk::Geometry3D*>(input->GetGeometry()->Clone().GetPointer()));
00108 
00109   output->SetPropertyList(input->GetPropertyList()->Clone());    
00110 
00111   m_TimeOfHeaderInitialization.Modified();
00112 }
00113 
00114 template < typename TPixel, unsigned int VImageDimension >
00115 void mitk::_InternalComputeClippedImage(itk::Image<TPixel, VImageDimension>* inputItkImage, mitk::GeometryClipImageFilter* geometryClipper, const mitk::Geometry2D* clippingGeometry2D)
00116 {
00117   typedef itk::Image<TPixel, VImageDimension> ItkInputImageType;
00118   typedef itk::Image<TPixel, VImageDimension> ItkOutputImageType;
00119 
00120   typedef itk::ImageRegionConstIteratorWithIndex< ItkInputImageType > ItkInputImageIteratorType;
00121   typedef itk::ImageRegionIteratorWithIndex< ItkOutputImageType > ItkOutputImageIteratorType;
00122 
00123   typename mitk::ImageToItk<ItkOutputImageType>::Pointer outputimagetoitk = mitk::ImageToItk<ItkOutputImageType>::New();
00124   outputimagetoitk->SetInput(geometryClipper->m_OutputTimeSelector->GetOutput());
00125   outputimagetoitk->Update();
00126   typename ItkOutputImageType::Pointer outputItkImage = outputimagetoitk->GetOutput();
00127 
00128   // create the iterators
00129   typename ItkInputImageType::RegionType inputRegionOfInterest = inputItkImage->GetLargestPossibleRegion();
00130   ItkInputImageIteratorType  inputIt( inputItkImage, inputRegionOfInterest );
00131   ItkOutputImageIteratorType outputIt( outputItkImage, inputRegionOfInterest );
00132 
00133   typename ItkOutputImageType::PixelType outsideValue;
00134   if(geometryClipper->m_AutoOutsideValue)
00135     outsideValue = itk::NumericTraits<typename ItkOutputImageType::PixelType>::min();
00136   else
00137     outsideValue = (typename ItkOutputImageType::PixelType) geometryClipper->m_OutsideValue;
00138 
00139   mitk::Geometry3D* inputGeometry = geometryClipper->m_InputTimeSelector->GetOutput()->GetGeometry();
00140   typedef itk::Index<VImageDimension> IndexType;
00141   Point3D indexPt; indexPt.Fill(0);
00142   int i, dim=IndexType::GetIndexDimension();
00143   Point3D pointInMM;
00144   bool above = geometryClipper->m_ClipPartAboveGeometry;
00145   bool labelBothSides = geometryClipper->GetLabelBothSides();
00146 
00147   if (geometryClipper->GetAutoOrientLabels())
00148   {
00149     Point3D leftMostPoint;
00150     leftMostPoint.Fill( std::numeric_limits<float>::min() / 2.0 );
00151     if(clippingGeometry2D->IsAbove(pointInMM) != above)
00152       {
00153       // invert meaning of above --> left is always the "above" side
00154       above = !above;
00155     MITK_INFO << leftMostPoint << " is BELOW geometry. Inverting meaning of above" << std::endl;
00156     }
00157     else
00158     MITK_INFO << leftMostPoint << " is above geometry" << std::endl;
00159   }
00160 
00161   typename ItkOutputImageType::PixelType aboveLabel = (typename ItkOutputImageType::PixelType)geometryClipper->GetAboveGeometryLabel();
00162   typename ItkOutputImageType::PixelType belowLabel = (typename ItkOutputImageType::PixelType)geometryClipper->GetBelowGeometryLabel();
00163 
00164   for ( inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
00165   {
00166     if((typename ItkOutputImageType::PixelType)inputIt.Get() == outsideValue)
00167     {
00168       outputIt.Set(outsideValue);
00169     }
00170     else
00171     {
00172       for(i=0;i<dim;++i)
00173         indexPt[i]=(mitk::ScalarType)inputIt.GetIndex()[i];
00174       inputGeometry->IndexToWorld(indexPt, pointInMM);
00175       if(clippingGeometry2D->IsAbove(pointInMM) == above)
00176       {
00177         if ( labelBothSides )
00178           outputIt.Set( aboveLabel );
00179         else
00180           outputIt.Set( outsideValue );
00181       }
00182       else
00183       {
00184         if ( labelBothSides)
00185           outputIt.Set( belowLabel );
00186         else
00187           outputIt.Set( inputIt.Get() );
00188       }
00189     }
00190   }
00191 }
00192 
00193 #include "mitkImageAccessByItk.h"
00194 
00195 void mitk::GeometryClipImageFilter::GenerateData()
00196 {
00197   Image::ConstPointer input = this->GetInput();
00198   Image::Pointer output = this->GetOutput();
00199 
00200   if((output->IsInitialized()==false) || (m_ClippingGeometry.IsNull()))
00201     return;
00202 
00203   const Geometry2D * clippingGeometryOfCurrentTimeStep = NULL;
00204 
00205   if(m_TimeSlicedClippingGeometry.IsNull())
00206   {
00207     clippingGeometryOfCurrentTimeStep = dynamic_cast<const Geometry2D*>(m_ClippingGeometry.GetPointer());
00208   }
00209   else
00210   {
00211     clippingGeometryOfCurrentTimeStep = dynamic_cast<const Geometry2D*>(m_TimeSlicedClippingGeometry->GetGeometry3D(0));
00212   }
00213 
00214   if(clippingGeometryOfCurrentTimeStep == NULL)
00215     return;
00216 
00217   m_InputTimeSelector->SetInput(input);
00218   m_OutputTimeSelector->SetInput(this->GetOutput());
00219 
00220   mitk::Image::RegionType outputRegion = output->GetRequestedRegion();
00221   const mitk::TimeSlicedGeometry *outputTimeGeometry = output->GetTimeSlicedGeometry();
00222   const mitk::TimeSlicedGeometry *inputTimeGeometry = input->GetTimeSlicedGeometry();
00223   ScalarType timeInMS;
00224 
00225   int timestep=0;
00226   int tstart=outputRegion.GetIndex(3);
00227   int tmax=tstart+outputRegion.GetSize(3);
00228 
00229   int t;
00230   for(t=tstart;t<tmax;++t)
00231   {
00232     timeInMS = outputTimeGeometry->TimeStepToMS( t );
00233 
00234     timestep = inputTimeGeometry->MSToTimeStep( timeInMS );
00235 
00236     m_InputTimeSelector->SetTimeNr(timestep);
00237     m_InputTimeSelector->UpdateLargestPossibleRegion();
00238     m_OutputTimeSelector->SetTimeNr(t);
00239     m_OutputTimeSelector->UpdateLargestPossibleRegion();
00240 
00241     if(m_TimeSlicedClippingGeometry.IsNotNull())
00242     {
00243       timestep = m_TimeSlicedClippingGeometry->MSToTimeStep( timeInMS );
00244       if(m_TimeSlicedClippingGeometry->IsValidTime(timestep) == false)
00245         continue;
00246 
00247       clippingGeometryOfCurrentTimeStep = dynamic_cast<const Geometry2D*>(m_TimeSlicedClippingGeometry->GetGeometry3D(timestep));
00248     }
00249 
00250     AccessByItk_2(m_InputTimeSelector->GetOutput(),_InternalComputeClippedImage,this,clippingGeometryOfCurrentTimeStep);
00251   }
00252 
00253   m_TimeOfHeaderInitialization.Modified();
00254 }