1 /*========================================================================= 2 * 3 * Copyright Insight Software Consortium 4 * 5 * Licensed under the Apache License, Version 2.0 (the "License"); 6 * you may not use this file except in compliance with the License. 7 * You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0.txt 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 * 17 *=========================================================================*/ 18 #ifndef itkLabelGeometryImageFilter_h 19 #define itkLabelGeometryImageFilter_h 20 21 #include "itkImageToImageFilter.h" 22 #include "itkNumericTraits.h" 23 #include "itkArray.h" 24 #include "itkSimpleDataObjectDecorator.h" 25 #include <map> 26 #include <vector> 27 #include "vnl/algo/vnl_symmetric_eigensystem.h" 28 #include "vnl/vnl_det.h" 29 #include "itkMath.h" 30 31 namespace itk 32 { 33 /** \class LabelGeometryImageFilter 34 * \brief Given a label map and an optional intensity image, compute 35 * geometric features. 36 * 37 * This filter enables the measurement of geometric features of all objects in 38 * a labeled ND volume. This labeled volume can represent, for instance, a 39 * medical image segmented into different anatomical structures or a microscope 40 * image segmented into individual cells. This filter is closely related to the 41 * itkLabelStatisticsImageFilter, which measures statistics of image regions 42 * defined by a labeled mask such as min, max, and mean intensity, intensity 43 * standard deviation, and bounding boxes. This filter, however, measures the 44 * geometry of the labeled regions themselves. 45 * 46 * It calculates features similar to the regionprops command of Matlab. The 47 * set of measurements that it enables include: volume, centroid, eigenvalues, 48 * eigenvectors, axes lenghts, eccentricity, elongation, orientation, bounding 49 * box, oriented bounding box, and rotation matrix. These features are based 50 * solely on the labeled mask itself. It also calculates integrated intensity 51 * and weighted centroid, which are measured on an intensity image under the 52 * labeled mask. These features represent the set of currently calculated 53 * features, but the framework of the filter is designed so that it can be 54 * easily expanded to measure a wide variety of other features. 55 * 56 * \authors Dirk Padfield and James Miller. 57 * 58 * This work is part of the National Alliance for Medical Image 59 * Computing (NAMIC), funded by the National Institutes of Health 60 * through the NIH Roadmap for Medical Research, Grant U54 EB005149. 61 * Information on the National Centers for Biomedical Computing 62 * can be obtained from http://commonfund.nih.gov/bioinformatics. 63 * 64 * This filter was contributed in the Insight Journal paper: 65 * "A Label Geometry Image Filter for Multiple Object Measurement" 66 * by Padfield D., Miller J 67 * http://www.insight-journal.org/browse/publication/301 68 * https://hdl.handle.net/1926/1493 69 * 70 * \ingroup ITKReview 71 * 72 * \wiki 73 * \wikiexample{ImageProcessing/LabelGeometryImageFilter,Get geometric properties of labeled regions in an image} 74 * \endwiki 75 */ 76 template< typename TLabelImage, typename TIntensityImage = TLabelImage > 77 class ITK_TEMPLATE_EXPORT LabelGeometryImageFilter: 78 public ImageToImageFilter< TLabelImage, TIntensityImage > 79 { 80 public: 81 ITK_DISALLOW_COPY_AND_ASSIGN(LabelGeometryImageFilter); 82 83 /** Standard Self type alias */ 84 using Self = LabelGeometryImageFilter; 85 using Superclass = ImageToImageFilter< TLabelImage, TIntensityImage >; 86 using Pointer = SmartPointer< Self >; 87 using ConstPointer = SmartPointer< const Self >; 88 89 /** Method for creation through the object factory. */ 90 itkNewMacro(Self); 91 92 /** Runtime information support. */ 93 itkTypeMacro(LabelGeometryImageFilter, ImageToImageFilter); 94 95 /** Image related type alias. */ 96 using IntensityImageType = TIntensityImage; 97 using InputImagePointer = typename TIntensityImage::Pointer; 98 using RegionType = typename TIntensityImage::RegionType; 99 using SizeType = typename TIntensityImage::SizeType; 100 using IndexType = typename TIntensityImage::IndexType; 101 using PixelType = typename TIntensityImage::PixelType; 102 103 /** Label image related type alias. */ 104 using LabelImageType = TLabelImage; 105 using LabelImagePointer = typename TLabelImage::Pointer; 106 using LabelRegionType = typename TLabelImage::RegionType; 107 using LabelSizeType = typename TLabelImage::SizeType; 108 using LabelIndexType = typename TLabelImage::IndexType; 109 using LabelPixelType = typename TLabelImage::PixelType; 110 using LabelPointType = typename TLabelImage::PointType; 111 112 /** Image related type alias. */ 113 static constexpr unsigned int ImageDimension = TLabelImage::ImageDimension; 114 115 /** Type to use for computations. */ 116 using RealType = typename NumericTraits< PixelType >::RealType; 117 118 /** Smart Pointer type to a DataObject. */ 119 using DataObjectPointer = typename DataObject::Pointer; 120 121 /** Type of DataObjects used for scalar outputs */ 122 using RealObjectType = SimpleDataObjectDecorator< RealType >; 123 124 /** Bounding Box-related type alias */ 125 using BoundingBoxType = itk::FixedArray< typename LabelIndexType::IndexValueType, 126 Self::ImageDimension *2 >; 127 using BoundingBoxFloatType = itk::FixedArray< float, 128 Self::ImageDimension *2 >; 129 130 // using BoundingBoxVerticesType = itk::FixedArray< 131 // LabelPointType,std::pow(2.0,Self::ImageDimension)>; 132 using BoundingBoxVerticesType = std::vector< LabelPointType >; 133 134 /** Axes Length-related type alias */ 135 using AxesLengthType = itk::FixedArray< RealType, Self::ImageDimension >; 136 137 /** Index array type alias */ 138 using IndexArrayType = itk::FixedArray< typename LabelIndexType::IndexValueType, 139 Self::ImageDimension >; 140 141 /** vector of labels */ 142 using LabelsType = std::vector< LabelPixelType >; 143 144 /** vector of indices */ 145 using LabelIndicesType = std::vector< LabelIndexType >; 146 147 /** Vector type */ 148 using VectorType = std::vector< double >; 149 150 /** Matrix type */ 151 using MatrixType = vnl_matrix< double >; 152 153 /** \class LabelGeometry 154 * \brief Geometry stored per label 155 * \ingroup ITKReview 156 */ 157 class LabelGeometry 158 { 159 public: 160 // default constructor LabelGeometry()161 LabelGeometry() 162 { 163 // initialized to the default values 164 this->m_Label = 0; 165 this->m_Sum = NumericTraits< RealType >::ZeroValue(); 166 167 const unsigned int imageDimension = Self::ImageDimension; 168 169 //m_BoundingBox.resize(imageDimension*2); 170 for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) 171 { 172 m_BoundingBox[i] = NumericTraits< typename IndexType::IndexValueType >::max(); 173 m_BoundingBox[i + 1] = NumericTraits< typename IndexType::IndexValueType >::NonpositiveMin(); 174 } 175 176 m_BoundingBoxVolume = 0; 177 m_BoundingBoxSize.Fill(0); 178 m_PixelIndices.clear(); 179 m_Centroid.Fill(0); 180 m_WeightedCentroid.Fill(0); 181 m_ZeroOrderMoment = 0; 182 m_FirstOrderRawMoments.Fill(0); 183 m_FirstOrderWeightedRawMoments.Fill(0); 184 m_Eigenvalues.resize(ImageDimension); 185 m_Eigenvalues.clear(); 186 m_Eigenvectors.set_size(ImageDimension, ImageDimension); 187 m_Eigenvectors.fill(0); 188 m_AxesLength.Fill(0); 189 m_Eccentricity = 1; 190 m_Elongation = 1; 191 m_Orientation = 0; 192 LabelPointType emptyPoint; 193 emptyPoint.Fill(0); 194 unsigned int numberOfVertices = 1 << ImageDimension; 195 m_OrientedBoundingBoxVertices.resize(numberOfVertices, emptyPoint); 196 m_OrientedBoundingBoxVolume = 0; 197 m_OrientedBoundingBoxSize.Fill(0); 198 m_OrientedLabelImage = LabelImageType::New(); 199 m_OrientedIntensityImage = IntensityImageType::New(); 200 m_OrientedBoundingBoxOrigin.Fill(0); 201 m_RotationMatrix.set_size(ImageDimension, ImageDimension); 202 m_RotationMatrix.fill(0.0); 203 204 m_SecondOrderRawMoments.set_size(ImageDimension, ImageDimension); 205 m_SecondOrderCentralMoments.set_size(ImageDimension, ImageDimension); 206 for ( unsigned int i = 0; i < ImageDimension; i++ ) 207 { 208 for ( unsigned int j = 0; j < ImageDimension; j++ ) 209 { 210 m_SecondOrderRawMoments(i, j) = 0; 211 m_SecondOrderCentralMoments(i, j) = 0; 212 } 213 } 214 } 215 216 LabelPixelType m_Label; 217 RealType m_Sum; 218 LabelPointType m_Centroid; 219 LabelPointType m_WeightedCentroid; 220 SizeValueType m_ZeroOrderMoment; 221 IndexArrayType m_FirstOrderRawMoments; 222 IndexArrayType m_FirstOrderWeightedRawMoments; 223 SizeValueType m_FirstOrderRawCrossMoment; 224 RealType m_FirstOrderCentralCrossMoment; 225 MatrixType m_SecondOrderRawMoments; 226 MatrixType m_SecondOrderCentralMoments; 227 VectorType m_Eigenvalues; 228 MatrixType m_Eigenvectors; 229 FixedArray< float, Self::ImageDimension > m_AxesLength; 230 RealType m_Eccentricity; 231 RealType m_Elongation; 232 RealType m_Orientation; 233 BoundingBoxType m_BoundingBox; 234 LabelSizeType m_BoundingBoxSize; 235 RealType m_BoundingBoxVolume; 236 LabelIndicesType m_PixelIndices; 237 BoundingBoxVerticesType m_OrientedBoundingBoxVertices; 238 RealType m_OrientedBoundingBoxVolume; 239 LabelPointType m_OrientedBoundingBoxSize; 240 typename LabelImageType::Pointer m_OrientedLabelImage; 241 typename IntensityImageType::Pointer m_OrientedIntensityImage; 242 MatrixType m_RotationMatrix; 243 LabelPointType m_OrientedBoundingBoxOrigin; 244 }; 245 246 /** Type of the map used to store data per label */ 247 // Map from the label to the class storing all of the geometry information. 248 using MapType = std::map< LabelPixelType, LabelGeometry >; 249 using MapIterator = typename std::map< LabelPixelType, LabelGeometry >::iterator; 250 using MapConstIterator = typename std::map< LabelPixelType, LabelGeometry >::const_iterator; 251 252 // Macros for enabling the calculation of additional features. 253 itkGetMacro(CalculatePixelIndices, bool); 254 itkBooleanMacro(CalculatePixelIndices); SetCalculatePixelIndices(const bool value)255 void SetCalculatePixelIndices(const bool value) 256 { 257 // CalculateOrientedBoundingBox, CalculateOrientedLabelImage, and 258 // CalculateOrientedIntensityImage all need CalculatePixelIndices to be 259 // turned 260 // on if they are turned on. So, CalculatePixelIndices cannot be 261 // turned off if any of these flags are turned on. 262 if ( value == false ) 263 { 264 if ( ( this->m_CalculateOrientedBoundingBox == true ) 265 || ( this->m_CalculateOrientedLabelRegions == true ) 266 || ( this->m_CalculateOrientedIntensityRegions == true ) ) 267 { 268 // We cannot change the value, so return. 269 return; 270 } 271 } 272 273 if ( this->m_CalculatePixelIndices != value ) 274 { 275 this->m_CalculatePixelIndices = value; 276 this->Modified(); 277 } 278 } 279 280 itkGetMacro(CalculateOrientedBoundingBox, bool); 281 itkBooleanMacro(CalculateOrientedBoundingBox); SetCalculateOrientedBoundingBox(const bool value)282 void SetCalculateOrientedBoundingBox(const bool value) 283 { 284 if ( this->m_CalculateOrientedBoundingBox != value ) 285 { 286 this->m_CalculateOrientedBoundingBox = value; 287 this->Modified(); 288 } 289 290 // CalculateOrientedBoundingBox needs 291 // CalculatePixelIndices to be turned on. 292 if ( value == true ) 293 { 294 this->SetCalculatePixelIndices(true); 295 } 296 } 297 298 itkGetMacro(CalculateOrientedLabelRegions, bool); 299 itkBooleanMacro(CalculateOrientedLabelRegions); SetCalculateOrientedLabelRegions(const bool value)300 void SetCalculateOrientedLabelRegions(const bool value) 301 { 302 if ( this->m_CalculateOrientedLabelRegions != value ) 303 { 304 this->m_CalculateOrientedLabelRegions = value; 305 this->Modified(); 306 307 // CalculateOrientedLabelImage needs 308 // CalculateOrientedBoundingBox to be turned on. 309 if ( value == true ) 310 { 311 SetCalculateOrientedBoundingBox(true); 312 } 313 } 314 } 315 316 itkGetMacro(CalculateOrientedIntensityRegions, bool); 317 itkBooleanMacro(CalculateOrientedIntensityRegions); SetCalculateOrientedIntensityRegions(const bool value)318 void SetCalculateOrientedIntensityRegions(const bool value) 319 { 320 if ( this->m_CalculateOrientedIntensityRegions != value ) 321 { 322 this->m_CalculateOrientedIntensityRegions = value; 323 this->Modified(); 324 325 // CalculateOrientedIntensityImage needs 326 // CalculateOrientedBoundingBox to be turned on. 327 if ( value == true ) 328 { 329 this->SetCalculateOrientedBoundingBox(true); 330 } 331 } 332 } 333 334 /** Set the intensity image */ SetIntensityInput(const TIntensityImage * input)335 void SetIntensityInput(const TIntensityImage *input) 336 { 337 // Process object is not const-correct so the const casting is required. 338 this->SetNthInput( 1, const_cast< TIntensityImage * >( input ) ); 339 } 340 341 /** Get the label image */ GetIntensityInput()342 const TIntensityImage * GetIntensityInput() const 343 { 344 return static_cast< TIntensityImage * >( const_cast< DataObject * >( this->ProcessObject::GetInput(1) ) ); 345 } 346 347 /** Does the specified label exist? Can only be called after 348 * a call to Update(). */ HasLabel(LabelPixelType label)349 bool HasLabel(LabelPixelType label) const 350 { 351 return m_LabelGeometryMapper.find(label) != m_LabelGeometryMapper.end(); 352 } 353 354 /** Get the number of labels used */ GetNumberOfObjects()355 SizeValueType GetNumberOfObjects() const 356 { 357 return m_LabelGeometryMapper.size(); 358 } 359 GetNumberOfLabels()360 SizeValueType GetNumberOfLabels() const 361 { 362 return this->GetNumberOfObjects(); 363 } 364 365 /** Get the labels that are in the image. */ GetLabels()366 std::vector< LabelPixelType > GetLabels() const 367 { 368 return m_AllLabels; 369 } 370 371 /** Return the all pixel indices for a label. */ 372 LabelIndicesType GetPixelIndices(LabelPixelType label) const; 373 374 /** Return the number of pixels for a label. This is the same as 375 * the volume and the zero order moment */ 376 SizeValueType GetVolume(LabelPixelType label) const; 377 378 /** Return the number of pixels for all labels. */ 379 //std::vector< SizeValueType > GetAllCounts() const; 380 381 /** Return the computed integrated pixel intensity for a label. */ 382 RealType GetIntegratedIntensity(LabelPixelType label) const; 383 384 /** Return the unweighted centroid for a label. */ 385 LabelPointType GetCentroid(LabelPixelType label) const; 386 387 /** Return the weighted centroid for a label. */ 388 LabelPointType GetWeightedCentroid(LabelPixelType label) const; 389 390 /** Return the eigenvalues as a vector. */ 391 VectorType GetEigenvalues(LabelPixelType label) const; 392 393 /** Return the eigenvectors as a matrix. */ 394 MatrixType GetEigenvectors(LabelPixelType label) const; 395 396 /** Return the axes length for a label. */ 397 AxesLengthType GetAxesLength(LabelPixelType label) const; 398 399 /** Return the minor axis length for a label. This is a convenience 400 * class that returns the shortest length from GetAxesLength. */ 401 RealType GetMinorAxisLength(LabelPixelType label) const; 402 403 /** Return the major axis length for a label. This is a convenience 404 * class that returns the longest length from GetAxesLength. */ 405 RealType GetMajorAxisLength(LabelPixelType label) const; 406 407 /** Return the eccentricity for a label. */ 408 RealType GetEccentricity(LabelPixelType label) const; 409 410 /** Return the elongation for a label. This is defined as the 411 * length of the major axis divided by the length of the minor axis. */ 412 RealType GetElongation(LabelPixelType label) const; 413 414 /** Return the orientation for a label defined in radians. */ 415 RealType GetOrientation(LabelPixelType label) const; 416 417 /** Return the computed bounding box for a label. 418 * This is organized in min/max pairs as [min(X), max(X), min(Y), 419 * max(Y), min(Z), max(Z),...] */ 420 BoundingBoxType GetBoundingBox(LabelPixelType label) const; 421 422 /** Return the volume of the bounding box. */ 423 RealType GetBoundingBoxVolume(LabelPixelType label) const; 424 425 /** Return the size of the bounding box. */ 426 LabelSizeType GetBoundingBoxSize(LabelPixelType label) const; 427 428 /** Return the oriented bounding box vertices. The order of the 429 * vertices corresponds with binary counting, where min is zero and 430 * max is one. For example, in 2D, binary counting will give 431 * [0,0],[0,1],[1,0],[1,1], which corresponds to 432 * [minX,minY],[minX,maxY],[maxX,minY],[maxX,maxY]. Each 433 * vertex is defined as an ND point. */ 434 BoundingBoxVerticesType GetOrientedBoundingBoxVertices(LabelPixelType label) const; 435 436 /** Return the volume of the oriented bounding box. */ 437 RealType GetOrientedBoundingBoxVolume(LabelPixelType label) const; 438 439 /** Return the size of the oriented bounding box. */ 440 LabelPointType GetOrientedBoundingBoxSize(LabelPixelType label) const; 441 442 /** Return the origin of the oriented bounding box. */ 443 LabelPointType GetOrientedBoundingBoxOrigin(LabelPixelType label) const; 444 445 /** Return the rotation matrix defined by the 446 * eigenvalues/eigenvectors. */ 447 MatrixType GetRotationMatrix(LabelPixelType label) const; 448 449 /** Return the region defined by the bounding box. */ 450 RegionType GetRegion(LabelPixelType label) const; 451 452 /** Return the label region defined by the oriented bounding box. */ 453 TLabelImage * GetOrientedLabelImage(LabelPixelType label) const; 454 455 /** Return the intensity region defined by the oriented bounding 456 * box. */ 457 TIntensityImage * GetOrientedIntensityImage(LabelPixelType label) const; 458 459 #ifdef ITK_USE_CONCEPT_CHECKING 460 // Begin concept checking 461 itkConceptMacro( InputHasNumericTraitsCheck, 462 ( Concept::HasNumericTraits< PixelType > ) ); 463 // End concept checking 464 #endif 465 466 protected: 467 LabelGeometryImageFilter(); 468 ~LabelGeometryImageFilter() override = default; 469 void PrintSelf(std::ostream & os, Indent indent) const override; 470 471 void GenerateData() override; 472 473 private: 474 bool CalculateOrientedBoundingBoxVertices(vnl_symmetric_eigensystem< double > eig, LabelGeometry & m_LabelGeometry); 475 476 bool m_CalculatePixelIndices; 477 bool m_CalculateOrientedBoundingBox; 478 bool m_CalculateOrientedLabelRegions; 479 bool m_CalculateOrientedIntensityRegions; 480 481 MapType m_LabelGeometryMapper; 482 LabelGeometry m_LabelGeometry; 483 LabelsType m_AllLabels; 484 }; // end of class 485 486 } // end namespace itk 487 488 #ifndef ITK_MANUAL_INSTANTIATION 489 #include "itkLabelGeometryImageFilter.hxx" 490 #endif 491 492 #endif 493