libinsighttoolkit3-dev 3.20.1-1 / usr / include / InsightToolkit / SpatialObject / itkImageMaskSpatialObject.txx

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    itkImageMaskSpatialObject.txx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkImageMaskSpatialObject_txx
#define __itkImageMaskSpatialObject_txx


#include "itkImageMaskSpatialObject.h"
#include "itkImageRegionConstIteratorWithIndex.h"

namespace itk
{

/** Constructor */
template< unsigned int TDimension>
ImageMaskSpatialObject< TDimension>
::ImageMaskSpatialObject()
{
  this->SetTypeName("ImageMaskSpatialObject");
  this->ComputeBoundingBox();
}

/** Destructor */
template< unsigned int TDimension>
ImageMaskSpatialObject< TDimension>
::~ImageMaskSpatialObject()
{
}


/** Test whether a point is inside or outside the object 
 *  For computational speed purposes, it is faster if the method does not
 *  check the name of the class and the current depth */ 
template< unsigned int TDimension >
bool 
ImageMaskSpatialObject< TDimension >
::IsInside( const PointType & point) const
{
  if( !this->GetBounds()->IsInside(point) )
    {
    return false;
    }

  if( !this->SetInternalInverseTransformToWorldToIndexTransform() )
    {
    return false;
    }

  PointType p = this->GetInternalInverseTransform()->TransformPoint(point);

  IndexType index;
  for(unsigned int i=0; i<TDimension; i++)
    {
    index[i] = static_cast<int>( p[i] );
    }

  const bool insideBuffer = this->GetImage()->GetBufferedRegion().IsInside( index );
     
  if( !insideBuffer )
    {
    return false;
    }

  const bool insideMask = (this->GetImage()->GetPixel(index) != NumericTraits<PixelType>::Zero);

  return insideMask;
}


/** Return true if the given point is inside the image */
template< unsigned int TDimension>
bool
ImageMaskSpatialObject< TDimension>
::IsInside( const PointType & point, unsigned int depth, char * name ) const
{
  if(name == NULL)
    {
    if(IsInside(point))
      {
      return true;
      }
    }
  else if(strstr(typeid(Self).name(), name))
    {
    if(IsInside(point))
      {
      return true;
      }
    }
  return SpatialObject<TDimension>::IsInside(point, depth, name);
}


template< unsigned int  TDimension >
typename ImageMaskSpatialObject< TDimension >::RegionType
ImageMaskSpatialObject< TDimension >
::GetAxisAlignedBoundingBoxRegion() const
{
  // We will use a slice iterator to iterate through slices orthogonal
  // to each of the axis of the image to find the bounding box. Each
  // slice iterator iterates from the outermost slice towards the image
  // center till it finds a mask pixel. For a 3D image, there will be six
  // slice iterators, iterating from the periphery inwards till the bounds
  // along each axes are found. The slice iterators save time and avoid 
  // having to walk the whole image. Since we are using slice iterators,
  // we will implement this only for 3D images.

  PixelType outsideValue = NumericTraits< PixelType >::Zero;
  RegionType region;
  
  ImagePointer image = this->GetImage();

  typedef typename SizeType::SizeValueType    SizeValueType;
  
  IndexType index;
  SizeType  size;

  if( ImageType::ImageDimension == 3)
    {
    
    for( unsigned int axis = 0; axis < ImageType::ImageDimension; axis++ )
      {
      // Two slice iterators along each axis...
      // Find the orthogonal planes for the slices
      unsigned int i, j;
      unsigned int direction[2];
      for (i = 0, j = 0; i < 3; ++i )
        {
        if (i != axis )
          {
          direction[j] = i;
          j++;
          }
        }
      

      // Create the forward iterator to find lower bound
      SliceIteratorType  fit(  image,  image->GetRequestedRegion() );
      fit.SetFirstDirection(  direction[1] );
      fit.SetSecondDirection( direction[0] );

      fit.GoToBegin();
      while( !fit.IsAtEnd() )
        {
        while( !fit.IsAtEndOfSlice() )
          {
          while( !fit.IsAtEndOfLine() )
            {
            if( fit.Get() != outsideValue )
              {
              index[axis] = fit.GetIndex()[axis];
              fit.GoToReverseBegin(); // skip to the end
              break;
              }
            ++fit;
            }
          fit.NextLine();
          }
        fit.NextSlice();
        }


      // Create the reverse iterator to find upper bound
      SliceIteratorType  rit(  image,  image->GetRequestedRegion() );
      rit.SetFirstDirection(  direction[1] );
      rit.SetSecondDirection( direction[0] );

      rit.GoToReverseBegin();
      while( !rit.IsAtReverseEnd() )
        {
        while( !rit.IsAtReverseEndOfSlice() )
          {
          while( !rit.IsAtReverseEndOfLine() )
            {
            if( rit.Get() != outsideValue )
              {
              size[axis] = rit.GetIndex()[axis] - index[axis] + 1;
              rit.GoToBegin(); //Skip to reverse end
              break;
              }
            --rit;
            }
          rit.PreviousLine();
          }
        rit.PreviousSlice();
        }
      }

    region.SetIndex( index );
    region.SetSize( size );
    }
  else
    {
    //itkExceptionMacro( << "ImageDimension must be 3!" );
    typedef ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
    IteratorType it( image, image->GetRequestedRegion() );
    it.GoToBegin();

    for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
      {
      index[ i ] = image->GetRequestedRegion().GetSize( i );
      size[ i ]  = image->GetRequestedRegion().GetIndex( i );
      }

    while( !it.IsAtEnd() )
      {
      if ( it.Get() != outsideValue )
        {
        IndexType tmpIndex = it.GetIndex();
        for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
          {

          if ( index[ i ] > tmpIndex[ i ] )
            {
            index[ i ] = tmpIndex[ i ];
            }

          const SizeValueType tmpSize = static_cast< SizeValueType >( tmpIndex[ i ] );

          if ( size[ i ]  < tmpSize )
            {
            size[ i ]  = tmpSize;
            }

          }
        }
      ++it;
      }

    for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
      {
      size[ i ] = size[ i ] - index[ i ] + 1;
      }
    region.SetIndex( index );
    region.SetSize( size );
  } // end else
  
  return region;
   
}


/** Print the object */
template< unsigned int TDimension >
void
ImageMaskSpatialObject< TDimension >
::PrintSelf( std::ostream& os, Indent indent ) const
{
  Superclass::PrintSelf(os,indent);
}


} // end namespace itk

#endif //__ImageMaskSpatialObject_txx