ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkLabelGeometryImageFilter.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkLabelGeometryImageFilter_h
00019 #define __itkLabelGeometryImageFilter_h
00020 
00021 #include "itkImageToImageFilter.h"
00022 #include "itkNumericTraits.h"
00023 #include "itkArray.h"
00024 #include "itkSimpleDataObjectDecorator.h"
00025 #include "itksys/hash_map.hxx"
00026 #include "itkFastMutexLock.h"
00027 #include <vector>
00028 #include "vnl/algo/vnl_symmetric_eigensystem.h"
00029 #include "vnl/vnl_det.h"
00030 #include "vnl/vnl_math.h"
00031 #include "vcl_cmath.h"
00032 
00033 namespace itk
00034 {
00078 template< class TLabelImage, class TIntensityImage = TLabelImage >
00079 class ITK_EXPORT LabelGeometryImageFilter:
00080   public ImageToImageFilter< TLabelImage, TIntensityImage >
00081 {
00082 public:
00084   typedef LabelGeometryImageFilter                           Self;
00085   typedef ImageToImageFilter< TLabelImage, TIntensityImage > Superclass;
00086   typedef SmartPointer< Self >                               Pointer;
00087   typedef SmartPointer< const Self >                         ConstPointer;
00088 
00090   itkNewMacro(Self);
00091 
00093   itkTypeMacro(LabelGeometryImageFilter, ImageToImageFilter);
00094 
00096   typedef TIntensityImage                      IntensityImageType;
00097   typedef typename TIntensityImage::Pointer    InputImagePointer;
00098   typedef typename TIntensityImage::RegionType RegionType;
00099   typedef typename TIntensityImage::SizeType   SizeType;
00100   typedef typename TIntensityImage::IndexType  IndexType;
00101   typedef typename TIntensityImage::PixelType  PixelType;
00102 
00104   typedef TLabelImage                      LabelImageType;
00105   typedef typename TLabelImage::Pointer    LabelImagePointer;
00106   typedef typename TLabelImage::RegionType LabelRegionType;
00107   typedef typename TLabelImage::SizeType   LabelSizeType;
00108   typedef typename TLabelImage::IndexType  LabelIndexType;
00109   typedef typename TLabelImage::PixelType  LabelPixelType;
00110   typedef typename TLabelImage::PointType  LabelPointType;
00111 
00113   itkStaticConstMacro(ImageDimension, unsigned int,
00114                       TLabelImage::ImageDimension);
00115 
00117   typedef typename NumericTraits< PixelType >::RealType RealType;
00118 
00120   typedef typename DataObject::Pointer DataObjectPointer;
00121 
00123   typedef SimpleDataObjectDecorator< RealType > RealObjectType;
00124 
00126   typedef itk::FixedArray< typename LabelIndexType::IndexValueType,
00127                            itkGetStaticConstMacro(ImageDimension) *2 > BoundingBoxType;
00128   typedef itk::FixedArray< float,
00129                            itkGetStaticConstMacro(ImageDimension) *2 > BoundingBoxFloatType;
00131 
00132   //typedef itk::FixedArray<
00133   // LabelPointType,vcl_pow(2.0,itkGetStaticConstMacro(ImageDimension))>
00134   // BoundingBoxVerticesType;
00135   typedef std::vector< LabelPointType > BoundingBoxVerticesType;
00136 
00138   typedef itk::FixedArray< RealType, itkGetStaticConstMacro(ImageDimension) > AxesLengthType;
00139 
00141   typedef itk::FixedArray< typename LabelIndexType::IndexValueType,
00142                            itkGetStaticConstMacro(ImageDimension) > IndexArrayType;
00143 
00145   typedef std::vector< LabelPixelType > LabelsType;
00146 
00148   typedef std::vector< LabelIndexType > LabelIndicesType;
00149 
00151   typedef std::vector< double > VectorType;
00152 
00154   typedef vnl_matrix< double > MatrixType;
00155 
00160   class LabelGeometry
00161   {
00162 public:
00163     // default constructor
00164     LabelGeometry()
00165     {
00166       // initialized to the default values
00167       this->m_Label = 0;
00168       this->m_Sum = NumericTraits< RealType >::Zero;
00169 
00170       const unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension);
00171 
00172       //m_BoundingBox.resize(imageDimension*2);
00173       for ( unsigned int i = 0; i < imageDimension * 2; i += 2 )
00174         {
00175         m_BoundingBox[i] = NumericTraits< typename IndexType::IndexValueType >::max();
00176         m_BoundingBox[i + 1] = NumericTraits< typename IndexType::IndexValueType >::NonpositiveMin();
00177         }
00178 
00179       m_BoundingBoxVolume = 0;
00180       m_BoundingBoxSize.Fill(0);
00181       m_PixelIndices.clear();
00182       m_Centroid.Fill(0);
00183       m_WeightedCentroid.Fill(0);
00184       m_ZeroOrderMoment = 0;
00185       m_FirstOrderRawMoments.Fill(0);
00186       m_FirstOrderWeightedRawMoments.Fill(0);
00187       m_Eigenvalues.resize(ImageDimension);
00188       m_Eigenvalues.clear();
00189       m_Eigenvectors.set_size(ImageDimension, ImageDimension);
00190       m_Eigenvectors.fill(0);
00191       m_AxesLength.Fill(0);
00192       m_Eccentricity = 1;
00193       m_Elongation = 1;
00194       m_Orientation = 0;
00195       LabelPointType emptyPoint;
00196       emptyPoint.Fill(0);
00197       unsigned int numberOfVertices = (unsigned int)vcl_pow( (double)2, (int)ImageDimension );
00198       m_OrientedBoundingBoxVertices.resize(numberOfVertices, emptyPoint);
00199       m_OrientedBoundingBoxVolume = 0;
00200       m_OrientedBoundingBoxSize.Fill(0);
00201       m_OrientedLabelImage = LabelImageType::New();
00202       m_OrientedIntensityImage = IntensityImageType::New();
00203       m_OrientedBoundingBoxOrigin.Fill(0);
00204       m_RotationMatrix.set_size(ImageDimension, ImageDimension);
00205       m_RotationMatrix.fill(0.0);
00206 
00207       m_SecondOrderRawMoments.set_size(ImageDimension, ImageDimension);
00208       m_SecondOrderCentralMoments.set_size(ImageDimension, ImageDimension);
00209       for ( unsigned int i = 0; i < ImageDimension; i++ )
00210         {
00211         for ( unsigned int j = 0; j < ImageDimension; j++ )
00212           {
00213           m_SecondOrderRawMoments(i, j) = 0;
00214           m_SecondOrderCentralMoments(i, j) = 0;
00215           }
00216         }
00217     }
00218 
00219     LabelPixelType                                              m_Label;
00220     RealType                                                    m_Sum;
00221     LabelPointType                                              m_Centroid;
00222     LabelPointType                                              m_WeightedCentroid;
00223     SizeValueType                                               m_ZeroOrderMoment;
00224     IndexArrayType                                              m_FirstOrderRawMoments;
00225     IndexArrayType                                              m_FirstOrderWeightedRawMoments;
00226     SizeValueType                                               m_FirstOrderRawCrossMoment;
00227     RealType                                                    m_FirstOrderCentralCrossMoment;
00228     MatrixType                                                  m_SecondOrderRawMoments;
00229     MatrixType                                                  m_SecondOrderCentralMoments;
00230     VectorType                                                  m_Eigenvalues;
00231     MatrixType                                                  m_Eigenvectors;
00232     FixedArray< float, itkGetStaticConstMacro(ImageDimension) > m_AxesLength;
00233     RealType                                                    m_Eccentricity;
00234     RealType                                                    m_Elongation;
00235     RealType                                                    m_Orientation;
00236     BoundingBoxType                                             m_BoundingBox;
00237     LabelSizeType                                               m_BoundingBoxSize;
00238     RealType                                                    m_BoundingBoxVolume;
00239     LabelIndicesType                                            m_PixelIndices;
00240     BoundingBoxVerticesType                                     m_OrientedBoundingBoxVertices;
00241     RealType                                                    m_OrientedBoundingBoxVolume;
00242     LabelPointType                                              m_OrientedBoundingBoxSize;
00243     typename LabelImageType::Pointer m_OrientedLabelImage;
00244     typename IntensityImageType::Pointer m_OrientedIntensityImage;
00245     MatrixType     m_RotationMatrix;
00246     LabelPointType m_OrientedBoundingBoxOrigin;
00247   };
00248 
00250   // Map from the label to the class storing all of the geometry information.
00251   typedef itksys::hash_map< LabelPixelType, LabelGeometry >                          MapType;
00252   typedef typename itksys::hash_map< LabelPixelType, LabelGeometry >::iterator       MapIterator;
00253   typedef typename itksys::hash_map< LabelPixelType, LabelGeometry >::const_iterator MapConstIterator;
00254 
00255   // Macros for enabling the calculation of additional features.
00256   itkGetMacro(CalculatePixelIndices, bool);
00257   itkBooleanMacro(CalculatePixelIndices);
00258   void SetCalculatePixelIndices(const bool value)
00259   {
00260     // CalculateOrientedBoundingBox, CalculateOrientedLabelImage, and
00261     // CalculateOrientedIntensityImage all need CalculatePixelIndices to be
00262     // turned
00263     // on if they are turned on.  So, CalculatePixelIndices cannot be
00264     // turned off if any of these flags are turned on.
00265     if ( value == false )
00266       {
00267       if ( ( this->m_CalculateOrientedBoundingBox == true )
00268            || ( this->m_CalculateOrientedLabelRegions == true )
00269            || ( this->m_CalculateOrientedIntensityRegions == true ) )
00270         {
00271         // We cannot change the value, so return.
00272         return;
00273         }
00274       }
00275 
00276     if ( this->m_CalculatePixelIndices != value )
00277       {
00278       this->m_CalculatePixelIndices = value;
00279       this->Modified();
00280       }
00281   }
00282 
00283   itkGetMacro(CalculateOrientedBoundingBox, bool);
00284   itkBooleanMacro(CalculateOrientedBoundingBox);
00285   void SetCalculateOrientedBoundingBox(const bool value)
00286   {
00287     if ( this->m_CalculateOrientedBoundingBox != value )
00288       {
00289       this->m_CalculateOrientedBoundingBox = value;
00290       this->Modified();
00291       }
00292 
00293     // CalculateOrientedBoundingBox needs
00294     // CalculatePixelIndices to be turned on.
00295     if ( value == true )
00296       {
00297       this->SetCalculatePixelIndices(true);
00298       }
00299   }
00300 
00301   itkGetMacro(CalculateOrientedLabelRegions, bool);
00302   itkBooleanMacro(CalculateOrientedLabelRegions);
00303   void SetCalculateOrientedLabelRegions(const bool value)
00304   {
00305     if ( this->m_CalculateOrientedLabelRegions != value )
00306       {
00307       this->m_CalculateOrientedLabelRegions = value;
00308       this->Modified();
00309 
00310       // CalculateOrientedLabelImage needs
00311       // CalculateOrientedBoundingBox to be turned on.
00312       if ( value == true )
00313         {
00314         SetCalculateOrientedBoundingBox(true);
00315         }
00316       }
00317   }
00318 
00319   itkGetMacro(CalculateOrientedIntensityRegions, bool);
00320   itkBooleanMacro(CalculateOrientedIntensityRegions);
00321   void SetCalculateOrientedIntensityRegions(const bool value)
00322   {
00323     if ( this->m_CalculateOrientedIntensityRegions != value )
00324       {
00325       this->m_CalculateOrientedIntensityRegions = value;
00326       this->Modified();
00327 
00328       // CalculateOrientedIntensityImage needs
00329       // CalculateOrientedBoundingBox to be turned on.
00330       if ( value == true )
00331         {
00332         this->SetCalculateOrientedBoundingBox(true);
00333         }
00334       }
00335   }
00336 
00338   void SetIntensityInput(const TIntensityImage *input)
00339   {
00340     // Process object is not const-correct so the const casting is required.
00341     this->SetNthInput( 1, const_cast< TIntensityImage * >( input ) );
00342   }
00343 
00345   const TIntensityImage * GetIntensityInput() const
00346   {
00347     return static_cast< TIntensityImage * >( const_cast< DataObject * >( this->ProcessObject::GetInput(1) ) );
00348   }
00349 
00352   bool HasLabel(LabelPixelType label) const
00353   {
00354     return m_LabelGeometryMapper.find(label) != m_LabelGeometryMapper.end();
00355   }
00356 
00358   SizeValueType GetNumberOfObjects() const
00359   {
00360     return m_LabelGeometryMapper.size();
00361   }
00362 
00363   SizeValueType GetNumberOfLabels() const
00364   {
00365     return this->GetNumberOfObjects();
00366   }
00367 
00369   std::vector< LabelPixelType > GetLabels() const
00370   {
00371     return m_AllLabels;
00372   }
00373 
00375   LabelIndicesType GetPixelIndices(LabelPixelType label) const;
00376 
00379   SizeValueType GetVolume(LabelPixelType label) const;
00380 
00382   //std::vector< SizeValueType > GetAllCounts() const;
00383 
00385   RealType GetIntegratedIntensity(LabelPixelType label) const;
00386 
00388   LabelPointType GetCentroid(LabelPixelType label) const;
00389 
00391   LabelPointType GetWeightedCentroid(LabelPixelType label) const;
00392 
00394   VectorType GetEigenvalues(LabelPixelType label) const;
00395 
00397   MatrixType GetEigenvectors(LabelPixelType label) const;
00398 
00400   AxesLengthType GetAxesLength(LabelPixelType label) const;
00401 
00404   RealType GetMinorAxisLength(LabelPixelType label) const;
00405 
00408   RealType GetMajorAxisLength(LabelPixelType label) const;
00409 
00411   RealType GetEccentricity(LabelPixelType label) const;
00412 
00415   RealType GetElongation(LabelPixelType label) const;
00416 
00418   RealType GetOrientation(LabelPixelType label) const;
00419 
00423   BoundingBoxType GetBoundingBox(LabelPixelType label) const;
00424 
00426   RealType GetBoundingBoxVolume(LabelPixelType label) const;
00427 
00429   LabelSizeType GetBoundingBoxSize(LabelPixelType label) const;
00430 
00437   BoundingBoxVerticesType GetOrientedBoundingBoxVertices(LabelPixelType label) const;
00438 
00440   RealType GetOrientedBoundingBoxVolume(LabelPixelType label) const;
00441 
00443   LabelPointType GetOrientedBoundingBoxSize(LabelPixelType label) const;
00444 
00446   LabelPointType GetOrientedBoundingBoxOrigin(LabelPixelType label) const;
00447 
00450   MatrixType GetRotationMatrix(LabelPixelType label) const;
00451 
00453   RegionType GetRegion(LabelPixelType label) const;
00454 
00456   TLabelImage * GetOrientedLabelImage(LabelPixelType label) const;
00457 
00460   TIntensityImage * GetOrientedIntensityImage(LabelPixelType label) const;
00461 
00462 #ifdef ITK_USE_CONCEPT_CHECKING
00463 
00464   itkConceptMacro( InputHasNumericTraitsCheck,
00465                    ( Concept::HasNumericTraits< PixelType > ) );
00466 
00468 #endif
00469 protected:
00470   LabelGeometryImageFilter();
00471   ~LabelGeometryImageFilter(){}
00472   void PrintSelf(std::ostream & os, Indent indent) const;
00474 
00475   void GenerateData();
00476 
00477 private:
00478   LabelGeometryImageFilter(const Self &); //purposely not implemented
00479   void operator=(const Self &);           //purposely not implemented
00480 
00481   bool CalculateOrientedBoundingBoxVertices(vnl_symmetric_eigensystem< double > eig, LabelGeometry & m_LabelGeometry);
00482 
00483   bool m_CalculatePixelIndices;
00484   bool m_CalculateOrientedBoundingBox;
00485   bool m_CalculateOrientedLabelRegions;
00486   bool m_CalculateOrientedIntensityRegions;
00487 
00488   MapType       m_LabelGeometryMapper;
00489   LabelGeometry m_LabelGeometry;
00490   LabelsType    m_AllLabels;
00491 
00492   SimpleFastMutexLock m_Mutex;
00493 }; // end of class
00494 
00495 template< class TLabelImage, class TIntensityImage >
00496 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::MatrixType CalculateRotationMatrix(
00497   vnl_symmetric_eigensystem< double > eig);
00498 
00499 template< class TLabelImage, class TIntensityImage, class TGenericImage >
00500 bool CalculateOrientedImage(
00501   LabelGeometryImageFilter< TLabelImage, TIntensityImage >  *filter,
00502   vnl_symmetric_eigensystem< double > eig,
00503   typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelGeometry & labelGeometry,
00504   bool useLabelImage = true);
00505 } // end namespace itk
00506 
00507 #ifndef ITK_MANUAL_INSTANTIATION
00508 #include "itkLabelGeometryImageFilter.hxx"
00509 #endif
00510 
00511 #endif
00512