itkLabelGeometryImageFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkLabelGeometryImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2010-01-24 14:09:29 $
  Version:   $Revision: 1.3 $

  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 __itkLabelGeometryImageFilter_h
#define __itkLabelGeometryImageFilter_h

#include "itkImageToImageFilter.h"
#include "itkNumericTraits.h"
#include "itkArray.h"
#include "itkSimpleDataObjectDecorator.h"
#include "itk_hash_map.h"
#include "itkFastMutexLock.h"
#include <vector>
#include <vnl/algo/vnl_symmetric_eigensystem.h>
#include <vnl/vnl_det.h>
#include <vnl/vnl_math.h>
#include <vcl_cmath.h>


namespace itk {

template<class TLabelImage, class TIntensityImage = TLabelImage>
class ITK_EXPORT LabelGeometryImageFilter : 
    public ImageToImageFilter<TLabelImage, TIntensityImage>
{
public:
  typedef LabelGeometryImageFilter                          Self;
  typedef ImageToImageFilter<TLabelImage,TIntensityImage>   Superclass;
  typedef SmartPointer<Self>                                Pointer;
  typedef SmartPointer<const Self>                          ConstPointer;

  itkNewMacro(Self);  

  itkTypeMacro(LabelGeometryImageFilter, ImageToImageFilter);

  typedef TIntensityImage                         IntensityImageType;
  typedef typename TIntensityImage::Pointer       InputImagePointer;
  typedef typename TIntensityImage::RegionType    RegionType;
  typedef typename TIntensityImage::SizeType      SizeType;
  typedef typename TIntensityImage::IndexType     IndexType;
  typedef typename TIntensityImage::PixelType     PixelType;

  typedef TLabelImage                             LabelImageType;
  typedef typename TLabelImage::Pointer           LabelImagePointer;
  typedef typename TLabelImage::RegionType        LabelRegionType;
  typedef typename TLabelImage::SizeType          LabelSizeType;
  typedef typename TLabelImage::IndexType         LabelIndexType;
  typedef typename TLabelImage::PixelType         LabelPixelType;
  typedef typename TLabelImage::PointType         LabelPointType;

  itkStaticConstMacro(ImageDimension, unsigned int,
                      TLabelImage::ImageDimension );

  typedef typename NumericTraits<PixelType>::RealType RealType;

  typedef typename DataObject::Pointer DataObjectPointer;

  typedef SimpleDataObjectDecorator<RealType> RealObjectType;

  typedef itk::FixedArray<typename LabelIndexType::IndexValueType,itkGetStaticConstMacro(ImageDimension)*2> BoundingBoxType;
  typedef itk::FixedArray< float,itkGetStaticConstMacro(ImageDimension)*2> BoundingBoxFloatType;

  //typedef itk::FixedArray< LabelPointType,vcl_pow(2.0,itkGetStaticConstMacro(ImageDimension))> BoundingBoxVerticesType;
  typedef std::vector< LabelPointType > BoundingBoxVerticesType;

  typedef itk::FixedArray<RealType,itkGetStaticConstMacro(ImageDimension)> AxesLengthType;

  typedef itk::FixedArray< typename LabelIndexType::IndexValueType, itkGetStaticConstMacro(ImageDimension) > IndexArrayType;

  typedef std::vector< LabelPixelType > LabelsType;

  typedef std::vector< LabelIndexType > LabelIndicesType;

  typedef std::vector< double > VectorType;

  typedef vnl_matrix<double> MatrixType;

  class LabelGeometry
    {
    public:
      // default constructor
      LabelGeometry()
      {
        // initialized to the default values
      this->m_Label = 0;
      this->m_Sum = NumericTraits<RealType>::Zero;

      const unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension);

      //m_BoundingBox.resize(imageDimension*2);
      for (unsigned int i = 0; i < imageDimension * 2; i += 2)
        {
        m_BoundingBox[i] = NumericTraits<ITK_TYPENAME IndexType::IndexValueType>::max();
        m_BoundingBox[i+1] = NumericTraits<ITK_TYPENAME IndexType::IndexValueType>::NonpositiveMin();
        }
        
      m_BoundingBoxVolume = 0;
      m_BoundingBoxSize.Fill(0);
      m_PixelIndices.clear();
      m_Centroid.Fill( 0 );
      m_WeightedCentroid.Fill( 0 );
      m_ZeroOrderMoment = 0;
      m_FirstOrderRawMoments.Fill( 0 );
      m_FirstOrderWeightedRawMoments.Fill( 0 );
      m_Eigenvalues.resize(ImageDimension);
      m_Eigenvalues.clear();
      m_Eigenvectors.set_size(ImageDimension,ImageDimension);
      m_Eigenvectors.fill(0);
      m_AxesLength.Fill( 0 );
      m_Eccentricity = 1;
      m_Elongation = 1;
      m_Orientation = 0;
      LabelPointType emptyPoint;
      emptyPoint.Fill( 0 );
      unsigned int numberOfVertices = (unsigned int)vcl_pow( (double)2, (int)ImageDimension );
      m_OrientedBoundingBoxVertices.resize(numberOfVertices,emptyPoint);
      m_OrientedBoundingBoxVolume = 0;
      m_OrientedBoundingBoxSize.Fill(0);
      m_OrientedLabelImage = LabelImageType::New();
      m_OrientedIntensityImage = IntensityImageType::New();
      m_OrientedBoundingBoxOrigin.Fill( 0 );
      m_RotationMatrix.set_size(ImageDimension,ImageDimension);
      m_RotationMatrix.fill( 0.0 );
      
      m_SecondOrderRawMoments.set_size(ImageDimension,ImageDimension);
      m_SecondOrderCentralMoments.set_size(ImageDimension,ImageDimension);
      for( unsigned int i = 0; i < ImageDimension; i++ )
        {
        for( unsigned int j = 0; j < ImageDimension; j++ )
          {
          m_SecondOrderRawMoments(i,j) = 0;
          m_SecondOrderCentralMoments(i,j) = 0;
          }
        }
      }

    LabelPixelType m_Label;
    RealType m_Sum;
    LabelPointType m_Centroid;
    LabelPointType m_WeightedCentroid;
    unsigned long m_ZeroOrderMoment;
    IndexArrayType m_FirstOrderRawMoments;
    IndexArrayType m_FirstOrderWeightedRawMoments;
    unsigned long  m_FirstOrderRawCrossMoment;
    RealType  m_FirstOrderCentralCrossMoment;
    MatrixType m_SecondOrderRawMoments;
    MatrixType m_SecondOrderCentralMoments;
    VectorType m_Eigenvalues;
    MatrixType m_Eigenvectors;
    FixedArray<float,itkGetStaticConstMacro(ImageDimension)> m_AxesLength;
    RealType m_Eccentricity;
    RealType m_Elongation;
    RealType m_Orientation;
    BoundingBoxType m_BoundingBox;
    LabelSizeType m_BoundingBoxSize;
    RealType m_BoundingBoxVolume;
    LabelIndicesType m_PixelIndices;
    BoundingBoxVerticesType m_OrientedBoundingBoxVertices;
    RealType m_OrientedBoundingBoxVolume;
    LabelPointType m_OrientedBoundingBoxSize;
    typename LabelImageType::Pointer m_OrientedLabelImage;
    typename IntensityImageType::Pointer m_OrientedIntensityImage;
    MatrixType m_RotationMatrix;
    LabelPointType m_OrientedBoundingBoxOrigin;
  };
  
  // Map from the label to the class storing all of the geometry information.
  typedef itk::hash_map<LabelPixelType, LabelGeometry> MapType;
  typedef typename itk::hash_map<LabelPixelType, LabelGeometry>::iterator MapIterator;
  typedef typename itk::hash_map<LabelPixelType, LabelGeometry>::const_iterator MapConstIterator;

  // Macros for enabling the calculation of additional features.
  itkGetMacro(CalculatePixelIndices, bool);
  itkBooleanMacro(CalculatePixelIndices);
  void SetCalculatePixelIndices( const bool value )
    {
    // CalculateOrientedBoundingBox, CalculateOrientedLabelImage, and
    // CalculateOrientedIntensityImage all need CalculatePixelIndices to be turned
    // on if they are turned on.  So, CalculatePixelIndices cannot be
    // turned off if any of these flags are turned on.
    if( value == false )
      {
      if( (this->m_CalculateOrientedBoundingBox == true) ||
          (this->m_CalculateOrientedLabelRegions == true) ||
          (this->m_CalculateOrientedIntensityRegions == true) )
        {
        // We cannot change the value, so return.
        return;
        }
      }
    
    if ( this->m_CalculatePixelIndices != value )
      {
      this->m_CalculatePixelIndices = value;
      this->Modified();
      }
    }
  
  itkGetMacro(CalculateOrientedBoundingBox, bool);
  itkBooleanMacro(CalculateOrientedBoundingBox);
  void SetCalculateOrientedBoundingBox( const bool value )
    {

    if (this->m_CalculateOrientedBoundingBox != value)
      {
      this->m_CalculateOrientedBoundingBox = value;
      this->Modified();
      }

    // CalculateOrientedBoundingBox needs
    // CalculatePixelIndices to be turned on.
    if( value == true )
      {
      this->SetCalculatePixelIndices( true );
      }

    }

  itkGetMacro(CalculateOrientedLabelRegions, bool);
  itkBooleanMacro(CalculateOrientedLabelRegions);
  void SetCalculateOrientedLabelRegions( const bool value )
    {
    if (this->m_CalculateOrientedLabelRegions != value)
      {
      this->m_CalculateOrientedLabelRegions = value;
      this->Modified();

      // CalculateOrientedLabelImage needs
      // CalculateOrientedBoundingBox to be turned on.
      if( value == true )
        {
        SetCalculateOrientedBoundingBox( true );
        }
      }
    }

  itkGetMacro(CalculateOrientedIntensityRegions, bool);
  itkBooleanMacro(CalculateOrientedIntensityRegions);
  void SetCalculateOrientedIntensityRegions( const bool value )
    {
    if (this->m_CalculateOrientedIntensityRegions != value)
      {
      this->m_CalculateOrientedIntensityRegions = value;
      this->Modified();

      // CalculateOrientedIntensityImage needs
      // CalculateOrientedBoundingBox to be turned on.
      if( value == true )
        {
        this->SetCalculateOrientedBoundingBox( true );
        }
      }
    }

  void SetIntensityInput(const TIntensityImage *input )
    {
    // Process object is not const-correct so the const casting is required.
    this->SetNthInput(1, const_cast<TIntensityImage *>(input) );
    }

  const TIntensityImage * GetIntensityInput() const
    {
    return static_cast<TIntensityImage*>(const_cast<DataObject *>(this->ProcessObject::GetInput(1)));
    }

  bool HasLabel(LabelPixelType label) const
    {
    return m_LabelGeometryMapper.find(label) != m_LabelGeometryMapper.end();
    }

  unsigned long GetNumberOfObjects() const
    {
    return m_LabelGeometryMapper.size();
    }
  unsigned long GetNumberOfLabels() const
    {
    return this->GetNumberOfObjects();
    }

  std::vector< LabelPixelType > GetLabels() const
    {
    return m_AllLabels;
    }

  LabelIndicesType GetPixelIndices( LabelPixelType label) const;

  unsigned long GetVolume(LabelPixelType label) const;

  //std::vector<unsigned long> GetAllCounts() const;

  RealType GetIntegratedIntensity(LabelPixelType label) const;

  LabelPointType GetCentroid( LabelPixelType label) const;

  LabelPointType GetWeightedCentroid( LabelPixelType label) const;

  VectorType GetEigenvalues (LabelPixelType label) const;

  MatrixType GetEigenvectors (LabelPixelType label) const;

  AxesLengthType GetAxesLength( LabelPixelType label) const;

  RealType GetMinorAxisLength( LabelPixelType label) const;

  RealType GetMajorAxisLength( LabelPixelType label) const;

  RealType GetEccentricity( LabelPixelType label) const;

  RealType GetElongation( LabelPixelType label) const;

  RealType GetOrientation( LabelPixelType label) const;

  BoundingBoxType GetBoundingBox(LabelPixelType label) const;

  RealType GetBoundingBoxVolume(LabelPixelType label) const;

  LabelSizeType GetBoundingBoxSize(LabelPixelType label) const;

  BoundingBoxVerticesType GetOrientedBoundingBoxVertices(LabelPixelType label) const;

  RealType GetOrientedBoundingBoxVolume(LabelPixelType label) const;

  LabelPointType GetOrientedBoundingBoxSize(LabelPixelType label) const;

  LabelPointType GetOrientedBoundingBoxOrigin(LabelPixelType label) const;

  MatrixType GetRotationMatrix(LabelPixelType label) const;

  RegionType GetRegion(LabelPixelType label) const;

  TLabelImage *GetOrientedLabelImage( LabelPixelType label) const;

  TIntensityImage * GetOrientedIntensityImage( LabelPixelType label) const;


#ifdef ITK_USE_CONCEPT_CHECKING

  itkConceptMacro(InputHasNumericTraitsCheck,
                  (Concept::HasNumericTraits<PixelType>));

#endif

protected:
  LabelGeometryImageFilter();
  ~LabelGeometryImageFilter(){};
  void PrintSelf(std::ostream& os, Indent indent) const;

  void GenerateData();

private:
  LabelGeometryImageFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

  bool CalculateOrientedBoundingBoxVertices(vnl_symmetric_eigensystem<double> eig, LabelGeometry & m_LabelGeometry);

  bool m_CalculatePixelIndices;
  bool m_CalculateOrientedBoundingBox;
  bool m_CalculateOrientedLabelRegions;
  bool m_CalculateOrientedIntensityRegions;

  MapType         m_LabelGeometryMapper;
  LabelGeometry   m_LabelGeometry;
  LabelsType      m_AllLabels;

  SimpleFastMutexLock m_Mutex;

}; // end of class

template<class TLabelImage, class TIntensityImage>
typename LabelGeometryImageFilter<TLabelImage,TIntensityImage>::MatrixType CalculateRotationMatrix(vnl_symmetric_eigensystem<double> eig);

template<class TLabelImage, class TIntensityImage, class TGenericImage>
bool CalculateOrientedImage(
  LabelGeometryImageFilter<TLabelImage, TIntensityImage>  *filter,
  vnl_symmetric_eigensystem<double> eig,
  typename LabelGeometryImageFilter<TLabelImage,TIntensityImage>::LabelGeometry & labelGeometry,
  bool useLabelImage = true);

} // end namespace itk
  
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkLabelGeometryImageFilter.txx"
#endif

#endif