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_hxx
19 #define itkLabelGeometryImageFilter_hxx
20
21 #include "itkLabelGeometryImageFilter.h"
22
23 #include "itkImageRegionConstIterator.h"
24 #include "itkImageRegionConstIteratorWithIndex.h"
25 #include "itkCastImageFilter.h"
26
27 #include "itkAffineTransform.h"
28 #include "itkResampleImageFilter.h"
29 #include "itkNearestNeighborInterpolateImageFunction.h"
30
31 namespace itk
32 {
33
34 namespace
35 {
36 //
37 // Private Helper functions
38 //
39 template< unsigned int NDimension >
40 vnl_matrix< double >
CalculateRotationMatrix(const vnl_symmetric_eigensystem<double> & eig)41 inline CalculateRotationMatrix(const vnl_symmetric_eigensystem< double > &eig)
42 {
43 vnl_matrix<double> rotationMatrix(NDimension, NDimension, 0);
44 for ( unsigned int i = 0; i < NDimension; i++ )
45 {
46 rotationMatrix.set_column( i, eig.get_eigenvector(i) );
47 }
48 // After vnl_symmetric_eigensystem, the columns of V are the eigenvectors,
49 // sorted by increasing eigenvalue, from most negative to most positive.
50 // First reorder the eigenvectors by decreasing eigenvalue.
51 rotationMatrix.fliplr();
52
53 // Next, check whether the determinant of the matrix is negative.
54 // If it is, then the vectors do not follow the right-hand rule. We
55 // can fix this by making one of them negative. Make the last
56 // eigenvector (with smallest eigenvalue) negative.
57 float matrixDet;
58 if ( NDimension == 2 )
59 {
60 matrixDet = vnl_det(rotationMatrix[0], rotationMatrix[1]);
61 }
62 else if ( NDimension == 3 )
63 {
64 matrixDet = vnl_det(rotationMatrix[0], rotationMatrix[1], rotationMatrix[2]);
65 }
66 else
67 {
68 matrixDet = 0.0f;
69 std::cerr << "ERROR: Determinant cannot be calculated for this dimension!" << std::endl;
70 }
71
72 if ( matrixDet < 0 )
73 {
74 rotationMatrix.set_column( NDimension - 1,
75 -rotationMatrix.get_column(NDimension - 1) );
76 }
77
78 // Transpose the matrix to yield the rotation matrix.
79 rotationMatrix.inplace_transpose();
80
81 return rotationMatrix;
82 }
83
84 template<typename TLabelImage, typename TIntensityImage, typename TInputImage >
85 inline bool
CalculateOrientedImage(const vnl_symmetric_eigensystem<double> & eig,typename LabelGeometryImageFilter<TLabelImage,TIntensityImage>::LabelGeometry & labelGeometry,bool useLabelImage,const TInputImage * inputImage)86 CalculateOrientedImage(
87 const vnl_symmetric_eigensystem< double > &eig,
88 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelGeometry & labelGeometry,
89 bool useLabelImage,
90 const TInputImage *inputImage)
91 {
92 static constexpr unsigned int Dimension = TInputImage::ImageDimension;
93
94 // CalculateOrientedBoundingBoxVertices needs to have already been
95 // called. This is taken care of by the flags.
96
97 // Rotate the original object to align with the coordinate
98 // system defined by the eigenvectors.
99 // Since the resampler moves from the output image to the input
100 // image, we take the transpose of the rotation matrix.
101 vnl_matrix<double> vnl_RotationMatrix = CalculateRotationMatrix< Dimension >(eig);
102 vnl_RotationMatrix.inplace_transpose();
103
104 // Set up the transform. Here the center of rotation is the
105 // centroid of the object, the rotation matrix is specified by the
106 // eigenvectors, and there is no translation.
107 using TransformType = itk::AffineTransform< double, Dimension >;
108 typename TransformType::Pointer transform = TransformType::New();
109 typename TransformType::MatrixType rotationMatrix(vnl_RotationMatrix);
110 typename TransformType::CenterType center;
111 typename TInputImage::PointType origin;
112 for( unsigned int i = 0; i < Dimension; i++ )
113 {
114 center[i] = labelGeometry.m_Centroid[i] * inputImage->GetSpacing()[i];
115 origin[i] = labelGeometry.m_OrientedBoundingBoxOrigin[i] * inputImage->GetSpacing()[i];
116 }
117 typename TransformType::OutputVectorType translation;
118 translation.Fill(0);
119 transform->SetCenter(center);
120 transform->SetTranslation(translation);
121 transform->SetMatrix(rotationMatrix);
122
123 using ResampleFilterType = itk::ResampleImageFilter< TInputImage, TIntensityImage >;
124 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
125
126 // The m_OrientedBoundingBoxSize is specified to float precision.
127 // Here we need an integer size large enough to contain all of the points, so
128 // we take the ceil of it.
129 // We also need to ensure that that bounding box is not outside of
130 // the image bounds.
131 typename ResampleFilterType::SizeType boundingBoxSize;
132 for ( unsigned int i = 0; i < Dimension; i++ )
133 {
134 boundingBoxSize[i] = ( typename ResampleFilterType::SizeType::SizeValueType )std::ceil(
135 labelGeometry.m_OrientedBoundingBoxSize[i]);
136 }
137
138 resampler->SetTransform(transform);
139 resampler->SetSize(boundingBoxSize);
140 resampler->SetOutputSpacing( inputImage->GetSpacing() );
141 resampler->SetOutputOrigin( origin );
142 resampler->SetInput(inputImage);
143
144 if ( useLabelImage )
145 {
146 // Set up the interpolator.
147 // Use a nearest neighbor interpolator since these are labeled images.
148 using InterpolatorType = itk::NearestNeighborInterpolateImageFunction< TInputImage, double >;
149 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
150 resampler->SetInterpolator(interpolator);
151 }
152 else
153 {
154 // Set up the interpolator.
155 // Use a linear interpolator since these are intensity images.
156 using InterpolatorType = itk::LinearInterpolateImageFunction< TInputImage, double >;
157 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
158 resampler->SetInterpolator(interpolator);
159 }
160 resampler->Update();
161 labelGeometry.m_OrientedIntensityImage->Graft( resampler->GetOutput() );
162
163 return true;
164 }
165
166 }
167
168 template< typename TLabelImage, typename TIntensityImage >
169 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
LabelGeometryImageFilter()170 ::LabelGeometryImageFilter()
171 {
172 this->SetNumberOfRequiredInputs(1);
173 m_CalculatePixelIndices = false;
174 m_CalculateOrientedBoundingBox = false;
175 m_CalculateOrientedLabelRegions = false;
176 m_CalculateOrientedIntensityRegions = false;
177 }
178
179 template< typename TLabelImage, typename TIntensityImage >
180 void
181 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GenerateData()182 ::GenerateData()
183 {
184 LabelPixelType label;
185
186 // Iterator over the label image.
187 ImageRegionConstIterator< TLabelImage > labelIt ( this->GetInput(),
188 this->GetInput()->GetBufferedRegion() );
189
190 using ImageIteratorWithIndexType = ImageRegionConstIteratorWithIndex< TLabelImage >;
191
192 // Iterator over the mapping from labels to geometry values.
193 MapIterator mapIt;
194
195 // begin with empty m_LabelGeometryMapper and m_AllLabels
196 m_LabelGeometryMapper.clear();
197 m_AllLabels.clear();
198
199 // Do the work
200 while ( !labelIt.IsAtEnd() )
201 {
202 label = labelIt.Get();
203
204 mapIt = m_LabelGeometryMapper.find(label);
205
206 // Is the label already in the mapper?
207 // If not, create a new geometry object.
208 if ( mapIt == m_LabelGeometryMapper.end() )
209 {
210 using MapValueType = typename MapType::value_type;
211
212 // map insert is defined as: pair<iterator, bool> insert(const value_type&
213 // x)
214 mapIt = m_LabelGeometryMapper.insert( MapValueType( label, LabelGeometry() ) ).first;
215 }
216
217 // Update the geometry values.
218
219 // LABEL
220 ( *mapIt ).second.m_Label = label;
221
222 // BOUNDING BOX
223 // The bounding box is defined in (min, max) pairs, such as
224 // (xmin,xmax,ymin,ymax,zmin,zmax).
225 typename ImageIteratorWithIndexType::IndexType index = labelIt.GetIndex();
226 for ( unsigned int i = 0; i < ( 2 * ImageDimension ); i += 2 )
227 {
228 // Update min
229 if ( ( *mapIt ).second.m_BoundingBox[i] > index[i / 2] )
230 {
231 ( *mapIt ).second.m_BoundingBox[i] = index[i / 2];
232 }
233 // Update max
234 if ( ( *mapIt ).second.m_BoundingBox[i + 1] < index[i / 2] )
235 {
236 ( *mapIt ).second.m_BoundingBox[i + 1] = index[i / 2];
237 }
238 }
239
240 // VOLUME
241 ( *mapIt ).second.m_ZeroOrderMoment++;
242
243 for ( unsigned int i = 0; i < ImageDimension; i++ )
244 {
245 // FIRST ORDER RAW MOMENTS
246 ( *mapIt ).second.m_FirstOrderRawMoments[i] += index[i];
247 }
248
249 // SECOND ORDER RAW MOMENTS
250 // Even for ND, the second order moments can be found from just
251 // two nested loops since second order moments consider only
252 // interactions between pairs of indices.
253 for ( unsigned int i = 0; i < ImageDimension; i++ )
254 {
255 // It is only necessary to fill in half of the matrix since it is
256 // symmetric.
257 for ( unsigned int j = 0; j < ImageDimension; j++ )
258 {
259 ( *mapIt ).second.m_SecondOrderRawMoments(i, j) += index[i] * index[j];
260 }
261 }
262
263 if ( m_CalculatePixelIndices == true )
264 {
265 // Pixel location list
266 ( *mapIt ).second.m_PixelIndices.push_back(index);
267 }
268
269 ++labelIt;
270 }
271
272 const TIntensityImage *intensityImage = this->GetIntensityInput();
273
274 // If an intensity image is defined, we can also calculate further
275 // features.
276 if ( intensityImage )
277 {
278 RealType value;
279
280 // Iterator over the intensity image.
281 using IntensityImageIteratorType = ImageRegionConstIteratorWithIndex< TIntensityImage >;
282
283 IntensityImageIteratorType it( intensityImage, intensityImage->GetBufferedRegion() );
284
285 typename IntensityImageIteratorType::IndexType index;
286
287 labelIt.GoToBegin();
288
289 while ( !it.IsAtEnd() )
290 {
291 label = labelIt.Get();
292 mapIt = m_LabelGeometryMapper.find(label);
293
294 value = static_cast< RealType >( it.Get() );
295 index = it.GetIndex();
296
297 // INTEGRATED PIXEL VALUE
298 ( *mapIt ).second.m_Sum += value;
299
300 for ( unsigned int i = 0; i < ImageDimension; i++ )
301 {
302 // FIRST ORDER WEIGHTED RAW MOMENTS
303 ( *mapIt ).second.m_FirstOrderWeightedRawMoments[i] += index[i]
304 * ( typename LabelIndexType::IndexValueType )value;
305 }
306
307 ++it;
308 ++labelIt;
309 }
310 }
311
312 // If there is no intensity input defined, the oriented
313 // intensity regions cannot be calculated.
314 if ( !intensityImage )
315 {
316 if ( m_CalculateOrientedIntensityRegions )
317 {
318 std::cerr
319 << "ERROR: An input intensity image must be used in order to calculate the oriented intensity image."
320 << std::endl;
321 }
322 m_CalculateOrientedIntensityRegions = false;
323 }
324
325 // We need to add to the second order moment the second order
326 // moment of a pixel. This can be derived analytically. The first
327 // order moment of a pixel can be shown to be 0, and the first order
328 // cross moment can also be shown to be 0. The second order moment
329 // can be shown to be 1/12.
330 float pixelSecondOrderCentralMoment = 1.0f / 12.0f;
331
332 // Now that the m_LabelGeometryMapper has been updated for all
333 // pixels in the image, we can calculate other geometrical values.
334 // Loop through all labels of the image.
335 for ( mapIt = m_LabelGeometryMapper.begin(); mapIt != m_LabelGeometryMapper.end(); mapIt++ )
336 {
337 // Update the bounding box measurements.
338 ( *mapIt ).second.m_BoundingBoxVolume = 1;
339 for ( unsigned int i = 0; i < ImageDimension; i++ )
340 {
341 ( *mapIt ).second.m_BoundingBoxSize[i] =
342 ( *mapIt ).second.m_BoundingBox[2 * i + 1] - ( *mapIt ).second.m_BoundingBox[2 * i] + 1;
343 ( *mapIt ).second.m_BoundingBoxVolume = ( *mapIt ).second.m_BoundingBoxVolume
344 * ( *mapIt ).second.m_BoundingBoxSize[i];
345 }
346
347 for ( unsigned int i = 0; i < ImageDimension; i++ )
348 {
349 // Normalize the centroid sum by the count to get the centroid.
350 ( *mapIt ).second.m_Centroid[i] =
351 static_cast< typename LabelPointType::ValueType >( ( *mapIt ).second.m_FirstOrderRawMoments[i] )
352 / ( *mapIt ).second.m_ZeroOrderMoment;
353
354 // This is the weighted sum. It only calculates correctly if
355 // the intensity image is defined.
356 if ( !intensityImage )
357 {
358 ( *mapIt ).second.m_WeightedCentroid[i] = 0.0;
359 }
360 else
361 {
362 ( *mapIt ).second.m_WeightedCentroid[i] =
363 static_cast< typename LabelPointType::ValueType >( ( *mapIt ).second.m_FirstOrderWeightedRawMoments[i] )
364 / ( *mapIt ).second.m_Sum;
365 }
366 }
367
368 // Using the raw moments, we can calculate the central moments.
369 MatrixType normalizedSecondOrderCentralMoments(ImageDimension, ImageDimension, 0);
370 for ( unsigned int i = 0; i < ImageDimension; i++ )
371 {
372 for ( unsigned int j = 0; j < ImageDimension; j++ )
373 {
374 normalizedSecondOrderCentralMoments(i,
375 j) =
376 ( ( *mapIt ).second.m_SecondOrderRawMoments(i,
377 j) ) / ( ( *mapIt ).second.m_ZeroOrderMoment )
378 - ( *mapIt ).second.m_Centroid[i]
379 * ( *mapIt ).second.m_Centroid[j];
380 // We need to add to the second order moment the second order
381 // moment of a pixel. This can be derived analytically.
382 if ( i == j )
383 {
384 normalizedSecondOrderCentralMoments(i, j) += pixelSecondOrderCentralMoment;
385 }
386 }
387 }
388
389 // Compute the eigenvalues/eigenvectors of the covariance matrix.
390 // The result is stored in increasing eigenvalues with
391 // corresponding eigenvectors.
392 vnl_symmetric_eigensystem< double > eig(normalizedSecondOrderCentralMoments);
393
394 // Calculate the eigenvalues/eigenvectors
395 VectorType eigenvalues(ImageDimension, 0);
396 MatrixType eigenvectors(ImageDimension, ImageDimension, 0);
397 for ( unsigned int i = 0; i < ImageDimension; i++ )
398 {
399 eigenvectors.set_column( i, eig.get_eigenvector(i) );
400 eigenvalues[i] = eig.get_eigenvalue(i);
401 }
402 ( *mapIt ).second.m_Eigenvalues = eigenvalues;
403 ( *mapIt ).second.m_Eigenvectors = eigenvectors;
404
405 itk::FixedArray< float, ImageDimension > axesLength;
406 for ( unsigned int i = 0; i < ImageDimension; i++ )
407 {
408 axesLength[i] = 4 * std::sqrt(eigenvalues[i]);
409 }
410 ( *mapIt ).second.m_AxesLength = axesLength;
411
412 // The following three features are currently only meaningful in 2D.
413 ( *mapIt ).second.m_Eccentricity = std::sqrt( ( eigenvalues[ImageDimension-1] - eigenvalues[0] ) / eigenvalues[ImageDimension-1] );
414 ( *mapIt ).second.m_Elongation = axesLength[ImageDimension-1] / axesLength[0];
415 RealType orientation = std::atan2(eig.get_eigenvector(ImageDimension-1)[1], eig.get_eigenvector(ImageDimension-1)[0]);
416 // Change the orientation from being between -pi to pi to being from 0 to pi.
417 // We can add pi because the orientation of the major axis is symmetric about the origin.
418 ( *mapIt ).second.m_Orientation = orientation < 0.0 ? orientation + itk::Math::pi : orientation;
419
420 if ( m_CalculateOrientedBoundingBox == true )
421 {
422 // Calculate the oriented bounding box using the eigenvectors.
423 CalculateOrientedBoundingBoxVertices(eig, ( *mapIt ).second);
424 }
425 if ( m_CalculateOrientedLabelRegions == true )
426 {
427 CalculateOrientedImage<TLabelImage, TIntensityImage>(
428 eig, ( *mapIt ).second, true, this->GetInput() );
429 }
430 if ( m_CalculateOrientedIntensityRegions == true )
431 {
432 // If there is no intensity input defined, the oriented
433 // intensity regions cannot be calculated.
434 if ( this->GetIntensityInput() )
435 {
436 CalculateOrientedImage<TLabelImage, TIntensityImage>(
437 eig, ( *mapIt ).second, false, this->GetIntensityInput() );
438 }
439 }
440
441 m_AllLabels.push_back( ( *mapIt ).first );
442 }
443 }
444
445 template< typename TLabelImage, typename TIntensityImage >
446 bool
447 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
CalculateOrientedBoundingBoxVertices(vnl_symmetric_eigensystem<double> eig,LabelGeometry & labelGeometry)448 ::CalculateOrientedBoundingBoxVertices(vnl_symmetric_eigensystem< double > eig, LabelGeometry & labelGeometry)
449 {
450 // Calculate the oriented bounding box using the eigenvectors.
451 // For each label, the pixels are rotated to the new coordinate
452 // system defined by the eigenvectors. The bounding boxes are
453 // calculated in this space, and then they are rotated back to the
454 // original coordinate system to define the oriented bounding boxes.
455 // The reverse rotation is the transpose of the rotation matrix.
456
457 // m_PixelIndices needs to have already been calculated. This is
458 // handled by the flags.
459
460 MatrixType rotationMatrix = CalculateRotationMatrix< TLabelImage::ImageDimension >(eig);
461 MatrixType inverseRotationMatrix = rotationMatrix.transpose();
462
463 labelGeometry.m_RotationMatrix = rotationMatrix;
464
465 // Convert the pixel locations to a vnl_matrix.
466 // Subtract the centroid of the region so that the rotation will
467 // be about the center of the region.
468 MatrixType pixelLocations(ImageDimension, labelGeometry.m_PixelIndices.size(), 0);
469 for ( unsigned int i = 0; i < labelGeometry.m_PixelIndices.size(); i++ )
470 {
471 for ( unsigned int j = 0; j < ImageDimension; j++ )
472 {
473 pixelLocations(j, i) = labelGeometry.m_PixelIndices[i][j] - labelGeometry.m_Centroid[j];
474 }
475 }
476
477 // Rotate the points by the rotation matrix from the eigenvectors.
478 MatrixType transformedPixelLocations = rotationMatrix * pixelLocations;
479
480 // Find the min and max of the point locations in this new
481 // coordinate system. These values will be float, so we use a
482 // BoundingBoxFloatType rather than BoundingBoxType.
483 // The bounding box order is [minX,maxX,minY,maxY,...]
484 BoundingBoxFloatType transformedBoundingBox;
485 // Initialize the bounding box values.
486 for ( unsigned int i = 0; i < ImageDimension * 2; i += 2 )
487 {
488 transformedBoundingBox[i] = NumericTraits< float >::max();
489 transformedBoundingBox[i + 1] = NumericTraits< float >::NonpositiveMin();
490 }
491
492 for ( unsigned int column = 0; column < transformedPixelLocations.columns(); column++ )
493 {
494 for ( unsigned int i = 0; i < ( 2 * ImageDimension ); i += 2 )
495 {
496 // Update min
497 if ( transformedBoundingBox[i] > transformedPixelLocations(i / 2, column) )
498 {
499 transformedBoundingBox[i] = transformedPixelLocations(i / 2, column);
500 }
501 // Update max
502 if ( transformedBoundingBox[i + 1] < transformedPixelLocations(i / 2, column) )
503 {
504 transformedBoundingBox[i + 1] = transformedPixelLocations(i / 2, column);
505 }
506 }
507 }
508
509 // Add 0.5 pixel buffers on each side of the bounding box to be sure to
510 // encompass the pixels and not cut through them.
511 for ( unsigned int i = 0; i < ( 2 * ImageDimension ); i += 2 )
512 {
513 // If the index corresponds with a min, subtract 0.5.
514 transformedBoundingBox[i] -= 0.5;
515 // Otherwise, add 0.5.
516 transformedBoundingBox[i + 1] += 0.5;
517 }
518
519 // The bounding box volume can be calculated directly from the
520 // transformed bounding box.
521 // Since 0.5 was already added to the border of the bounding box,
522 // it is not necessary to add 1 to the range to get the correct range.
523 labelGeometry.m_OrientedBoundingBoxVolume = 1;
524 for ( unsigned int i = 0; i < ( 2 * ImageDimension ); i += 2 )
525 {
526 labelGeometry.m_OrientedBoundingBoxSize[i / 2] = transformedBoundingBox[i + 1] - transformedBoundingBox[i];
527 labelGeometry.m_OrientedBoundingBoxVolume *= transformedBoundingBox[i + 1] - transformedBoundingBox[i];
528 }
529
530 // The bounding box cannot be transformed directly because an
531 // oriented bounding box cannot be specified simply by the min and
532 // max in each dimension. Instead, the oriented bounding box is
533 // specified by the points corresponding to the vertices of the
534 // oriented bounding box. We find these from the oriented
535 // bounding box and transform them back to the original coordinate frame.
536 // The order of the vertices corresponds with binary counting,
537 // where min is zero and max is one. For example, in 2D, binary
538 // counting will give [0,0],[0,1],[1,0],[1,1], which corresponds
539 // to [minX,minY],[minX,maxY],[maxX,minY],[maxX,maxY].
540 // Loop through each dimension of the bounding box and find all of the
541 // vertices.
542 unsigned int numberOfVertices = 1 << ImageDimension;
543 MatrixType transformedBoundingBoxVertices(ImageDimension, numberOfVertices, 0);
544 int val;
545 LabelIndexType binaryIndex;
546 int arrayIndex;
547 for ( unsigned int i = 0; i < numberOfVertices; i++ )
548 {
549 val = i;
550 for ( unsigned int j = 0; j < ImageDimension; j++ )
551 {
552 // This is the binary index as described above.
553 binaryIndex[j] = val % 2;
554 val = val / 2;
555 // This is the index into the transformedBoundingBox array
556 // corresponding to the binaryIndex.
557 arrayIndex = binaryIndex[j] + 2 * j;
558 transformedBoundingBoxVertices(j, i) = transformedBoundingBox[arrayIndex];
559 }
560 }
561
562 // Transform the transformed bounding box vertices back to the
563 // original coordinate system.
564 MatrixType orientedBoundingBoxVertices = inverseRotationMatrix * transformedBoundingBoxVertices;
565
566 // Add the centroid back to each of the vertices since it was
567 // subtracted when the points were rotated.
568 for ( unsigned int i = 0; i < orientedBoundingBoxVertices.columns(); i++ )
569 {
570 for ( unsigned int j = 0; j < ImageDimension; j++ )
571 {
572 orientedBoundingBoxVertices(j, i) += labelGeometry.m_Centroid[j];
573 // Copy the oriented bounding box vertices back to a vector of
574 // points for the mapper.
575 labelGeometry.m_OrientedBoundingBoxVertices[i][j] = orientedBoundingBoxVertices(j, i);
576 }
577 }
578
579 // Find the origin of the oriented bounding box.
580 for ( unsigned int i = 0; i < ImageDimension; i++ )
581 {
582 labelGeometry.m_OrientedBoundingBoxOrigin[i] = transformedBoundingBox[2 * i] + labelGeometry.m_Centroid[i];
583 }
584
585 return true;
586 }
587
588 template< typename TLabelImage, typename TIntensityImage >
589 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelIndicesType
590 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetPixelIndices(LabelPixelType label) const591 ::GetPixelIndices(LabelPixelType label) const
592 {
593 MapConstIterator mapIt;
594
595 mapIt = m_LabelGeometryMapper.find(label);
596 if ( mapIt == m_LabelGeometryMapper.end() )
597 {
598 // label does not exist, return a default value
599 LabelIndicesType emptyVector;
600 emptyVector.clear();
601 return emptyVector;
602 }
603 else
604 {
605 return ( *mapIt ).second.m_PixelIndices;
606 }
607 }
608
609 template< typename TLabelImage, typename TIntensityImage >
610 SizeValueType
611 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetVolume(LabelPixelType label) const612 ::GetVolume(LabelPixelType label) const
613 {
614 MapConstIterator mapIt;
615
616 mapIt = m_LabelGeometryMapper.find(label);
617 if ( mapIt == m_LabelGeometryMapper.end() )
618 {
619 // label does not exist, return a default value
620 return 0;
621 }
622 else
623 {
624 return ( *mapIt ).second.m_ZeroOrderMoment;
625 }
626 }
627
628 template< typename TLabelImage, typename TIntensityImage >
629 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
630 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetIntegratedIntensity(LabelPixelType label) const631 ::GetIntegratedIntensity(LabelPixelType label) const
632 {
633 MapConstIterator mapIt;
634
635 mapIt = m_LabelGeometryMapper.find(label);
636 if ( mapIt == m_LabelGeometryMapper.end() )
637 {
638 // label does not exist, return a default value
639 return NumericTraits< RealType >::ZeroValue();
640 }
641 else
642 {
643 return ( *mapIt ).second.m_Sum;
644 }
645 }
646
647 template< typename TLabelImage, typename TIntensityImage >
648 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelPointType
649 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetCentroid(LabelPixelType label) const650 ::GetCentroid(LabelPixelType label) const
651 {
652 MapConstIterator mapIt;
653
654 mapIt = m_LabelGeometryMapper.find(label);
655 if ( mapIt == m_LabelGeometryMapper.end() )
656 {
657 // label does not exist, return a default value
658 LabelPointType emptyCentroid;
659 emptyCentroid.Fill(NumericTraits< typename LabelPointType::ValueType >::ZeroValue());
660 return emptyCentroid;
661 }
662 else
663 {
664 return ( *mapIt ).second.m_Centroid;
665 }
666 }
667
668 template< typename TLabelImage, typename TIntensityImage >
669 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelPointType
670 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetWeightedCentroid(LabelPixelType label) const671 ::GetWeightedCentroid(LabelPixelType label) const
672 {
673 MapConstIterator mapIt;
674
675 mapIt = m_LabelGeometryMapper.find(label);
676 if ( mapIt == m_LabelGeometryMapper.end() )
677 {
678 // label does not exist, return a default value
679 LabelPointType emptyCentroid;
680 emptyCentroid.Fill(NumericTraits< typename LabelPointType::ValueType >::ZeroValue());
681 return emptyCentroid;
682 }
683 else
684 {
685 return ( *mapIt ).second.m_WeightedCentroid;
686 }
687 }
688
689 template< typename TLabelImage, typename TIntensityImage >
690 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::VectorType
691 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetEigenvalues(LabelPixelType label) const692 ::GetEigenvalues(LabelPixelType label) const
693 {
694 MapConstIterator mapIt;
695
696 mapIt = m_LabelGeometryMapper.find(label);
697 if ( mapIt == m_LabelGeometryMapper.end() )
698 {
699 // label does not exist, return a default value
700 VectorType emptyVector(ImageDimension, 0);
701 return emptyVector;
702 }
703 else
704 {
705 return ( *mapIt ).second.m_Eigenvalues;
706 }
707 }
708
709 template< typename TLabelImage, typename TIntensityImage >
710 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::MatrixType
711 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetEigenvectors(LabelPixelType label) const712 ::GetEigenvectors(LabelPixelType label) const
713 {
714 MapConstIterator mapIt;
715
716 mapIt = m_LabelGeometryMapper.find(label);
717 if ( mapIt == m_LabelGeometryMapper.end() )
718 {
719 // label does not exist, return a default value
720 MatrixType emptyMatrix(ImageDimension, ImageDimension, 0);
721 return emptyMatrix;
722 }
723 else
724 {
725 return ( *mapIt ).second.m_Eigenvectors;
726 }
727 }
728
729 template< typename TLabelImage, typename TIntensityImage >
730 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::AxesLengthType
731 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetAxesLength(LabelPixelType label) const732 ::GetAxesLength(LabelPixelType label) const
733 {
734 MapConstIterator mapIt;
735
736 mapIt = m_LabelGeometryMapper.find(label);
737 if ( mapIt == m_LabelGeometryMapper.end() )
738 {
739 // label does not exist, return a default value
740 LabelPointType emptyAxesLength;
741 emptyAxesLength.Fill(NumericTraits< typename AxesLengthType::ValueType >::ZeroValue());
742 return emptyAxesLength;
743 }
744 else
745 {
746 return ( *mapIt ).second.m_AxesLength;
747 }
748 }
749
750 template< typename TLabelImage, typename TIntensityImage >
751 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
752 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetMinorAxisLength(LabelPixelType label) const753 ::GetMinorAxisLength(LabelPixelType label) const
754 {
755 AxesLengthType axisLength = GetAxesLength(label);
756
757 return axisLength[0];
758 }
759
760 template< typename TLabelImage, typename TIntensityImage >
761 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
762 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetMajorAxisLength(LabelPixelType label) const763 ::GetMajorAxisLength(LabelPixelType label) const
764 {
765 AxesLengthType axisLength = GetAxesLength(label);
766
767 return axisLength[ImageDimension - 1];
768 }
769
770 template< typename TLabelImage, typename TIntensityImage >
771 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
772 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetEccentricity(LabelPixelType label) const773 ::GetEccentricity(LabelPixelType label) const
774 {
775 MapConstIterator mapIt;
776
777 mapIt = m_LabelGeometryMapper.find(label);
778 if ( mapIt == m_LabelGeometryMapper.end() )
779 {
780 // label does not exist, return a default value
781 return NumericTraits< RealType >::ZeroValue();
782 }
783 else
784 {
785 return ( *mapIt ).second.m_Eccentricity;
786 }
787 }
788
789 template< typename TLabelImage, typename TIntensityImage >
790 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
791 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetElongation(LabelPixelType label) const792 ::GetElongation(LabelPixelType label) const
793 {
794 MapConstIterator mapIt;
795
796 mapIt = m_LabelGeometryMapper.find(label);
797 if ( mapIt == m_LabelGeometryMapper.end() )
798 {
799 // label does not exist, return a default value
800 return NumericTraits< RealType >::ZeroValue();
801 }
802 else
803 {
804 return ( *mapIt ).second.m_Elongation;
805 }
806 }
807
808 template< typename TLabelImage, typename TIntensityImage >
809 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
810 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetOrientation(LabelPixelType label) const811 ::GetOrientation(LabelPixelType label) const
812 {
813 MapConstIterator mapIt;
814
815 mapIt = m_LabelGeometryMapper.find(label);
816 if ( mapIt == m_LabelGeometryMapper.end() )
817 {
818 // label does not exist, return a default value
819 return NumericTraits< RealType >::ZeroValue();
820 }
821 else
822 {
823 return ( *mapIt ).second.m_Orientation;
824 }
825 }
826
827 template< typename TLabelImage, typename TIntensityImage >
828 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::BoundingBoxType
829 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetBoundingBox(LabelPixelType label) const830 ::GetBoundingBox(LabelPixelType label) const
831 {
832 MapConstIterator mapIt;
833
834 mapIt = m_LabelGeometryMapper.find(label);
835 if ( mapIt == m_LabelGeometryMapper.end() )
836 {
837 BoundingBoxType emptyBox;
838 emptyBox.Fill(NumericTraits< typename BoundingBoxType::ValueType >::ZeroValue());
839 // label does not exist, return a default value
840 return emptyBox;
841 }
842 else
843 {
844 return ( *mapIt ).second.m_BoundingBox;
845 }
846 }
847
848 template< typename TLabelImage, typename TIntensityImage >
849 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
850 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetBoundingBoxVolume(LabelPixelType label) const851 ::GetBoundingBoxVolume(LabelPixelType label) const
852 {
853 MapConstIterator mapIt;
854
855 mapIt = m_LabelGeometryMapper.find(label);
856 if ( mapIt == m_LabelGeometryMapper.end() )
857 {
858 // label does not exist, return a default value
859 return NumericTraits< RealType >::ZeroValue();
860 }
861 else
862 {
863 return ( *mapIt ).second.m_BoundingBoxVolume;
864 }
865 }
866
867 template< typename TLabelImage, typename TIntensityImage >
868 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelSizeType
869 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetBoundingBoxSize(LabelPixelType label) const870 ::GetBoundingBoxSize(LabelPixelType label) const
871 {
872 MapConstIterator mapIt;
873
874 mapIt = m_LabelGeometryMapper.find(label);
875 if ( mapIt == m_LabelGeometryMapper.end() )
876 {
877 // label does not exist, return a default value
878 LabelSizeType emptySize;
879 emptySize.Fill(NumericTraits< typename LabelSizeType::SizeValueType >::ZeroValue());
880 return emptySize;
881 }
882 else
883 {
884 return ( *mapIt ).second.m_BoundingBoxSize;
885 }
886 }
887
888 template< typename TLabelImage, typename TIntensityImage >
889 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::BoundingBoxVerticesType
890 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetOrientedBoundingBoxVertices(LabelPixelType label) const891 ::GetOrientedBoundingBoxVertices(LabelPixelType label) const
892 {
893 unsigned int numberOfVertices = 1 << ImageDimension;
894 MapConstIterator mapIt = m_LabelGeometryMapper.find(label);
895 if ( mapIt == m_LabelGeometryMapper.end() )
896 {
897 // label does not exist, return a default value
898 LabelPointType emptyPoint;
899 emptyPoint.Fill(0);
900 BoundingBoxVerticesType emptyVertices;
901 emptyVertices.resize(numberOfVertices, emptyPoint);
902 return emptyVertices;
903 }
904 else
905 {
906 return ( *mapIt ).second.m_OrientedBoundingBoxVertices;
907 }
908 }
909
910 template< typename TLabelImage, typename TIntensityImage >
911 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RealType
912 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetOrientedBoundingBoxVolume(LabelPixelType label) const913 ::GetOrientedBoundingBoxVolume(LabelPixelType label) const
914 {
915 MapConstIterator mapIt;
916
917 mapIt = m_LabelGeometryMapper.find(label);
918 if ( mapIt == m_LabelGeometryMapper.end() )
919 {
920 // label does not exist, return a default value
921 return NumericTraits< RealType >::ZeroValue();
922 }
923 else
924 {
925 return ( *mapIt ).second.m_OrientedBoundingBoxVolume;
926 }
927 }
928
929 template< typename TLabelImage, typename TIntensityImage >
930 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelPointType
931 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetOrientedBoundingBoxSize(LabelPixelType label) const932 ::GetOrientedBoundingBoxSize(LabelPixelType label) const
933 {
934 MapConstIterator mapIt;
935
936 mapIt = m_LabelGeometryMapper.find(label);
937 if ( mapIt == m_LabelGeometryMapper.end() )
938 {
939 // label does not exist, return a default value
940 // LabelSizeType emptySize;
941 // emptySize.Fill( NumericTraits<LabelSizeType::SizeValueType>::ZeroValue());
942 // return emptySize;
943 LabelPointType emptySize;
944 emptySize.Fill(NumericTraits< typename LabelPointType::ValueType >::ZeroValue());
945 return emptySize;
946 }
947 else
948 {
949 return ( *mapIt ).second.m_OrientedBoundingBoxSize;
950 }
951 }
952
953 template< typename TLabelImage, typename TIntensityImage >
954 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::LabelPointType
955 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetOrientedBoundingBoxOrigin(LabelPixelType label) const956 ::GetOrientedBoundingBoxOrigin(LabelPixelType label) const
957 {
958 MapConstIterator mapIt;
959
960 mapIt = m_LabelGeometryMapper.find(label);
961 if ( mapIt == m_LabelGeometryMapper.end() )
962 {
963 // label does not exist, return a default value
964 LabelPointType emptySize;
965 emptySize.Fill(NumericTraits< typename LabelPointType::ValueType >::ZeroValue());
966 return emptySize;
967 }
968 else
969 {
970 return ( *mapIt ).second.m_OrientedBoundingBoxOrigin;
971 }
972 }
973
974 template< typename TLabelImage, typename TIntensityImage >
975 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::MatrixType
976 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetRotationMatrix(LabelPixelType label) const977 ::GetRotationMatrix(LabelPixelType label) const
978 {
979 MapConstIterator mapIt;
980
981 mapIt = m_LabelGeometryMapper.find(label);
982 if ( mapIt == m_LabelGeometryMapper.end() )
983 {
984 // label does not exist, return a default value
985 MatrixType emptyMatrix(ImageDimension, ImageDimension, 0);
986 return emptyMatrix;
987 }
988 else
989 {
990 return ( *mapIt ).second.m_RotationMatrix;
991 }
992 }
993
994 template< typename TLabelImage, typename TIntensityImage >
995 typename LabelGeometryImageFilter< TLabelImage, TIntensityImage >::RegionType
996 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetRegion(LabelPixelType label) const997 ::GetRegion(LabelPixelType label) const
998 {
999 MapConstIterator mapIt;
1000
1001 mapIt = m_LabelGeometryMapper.find(label);
1002
1003 if ( mapIt == m_LabelGeometryMapper.end() )
1004 {
1005 RegionType emptyRegion;
1006 // label does not exist, return a default value
1007 return emptyRegion;
1008 }
1009 else
1010 {
1011 BoundingBoxType bbox = this->GetBoundingBox(label);
1012 IndexType index;
1013 SizeType size;
1014
1015 unsigned int dimension = bbox.Size() / 2;
1016
1017 for ( unsigned int i = 0; i < dimension; i++ )
1018 {
1019 index[i] = bbox[2 * i];
1020 size[i] = bbox[2 * i + 1] - bbox[2 * i] + 1;
1021 }
1022 RegionType region;
1023 region.SetSize(size);
1024 region.SetIndex(index);
1025
1026 return region;
1027 }
1028 }
1029
1030 template< typename TLabelImage, typename TIntensityImage >
1031 TLabelImage *
1032 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
1033 ::GetOrientedLabelImage(LabelPixelType label) const
1034 {
1035 MapConstIterator mapIt;
1036
1037 mapIt = m_LabelGeometryMapper.find(label);
1038 if ( mapIt == m_LabelGeometryMapper.end() )
1039 {
1040 // label does not exist, return a default value
1041 return nullptr;
1042 }
1043 else
1044 {
1045 return ( *mapIt ).second.m_OrientedLabelImage;
1046 }
1047 }
1048
1049 template< typename TLabelImage, typename TIntensityImage >
1050 TIntensityImage *
1051 LabelGeometryImageFilter< TLabelImage, TIntensityImage >
GetOrientedIntensityImage(LabelPixelType label) const1052 ::GetOrientedIntensityImage(LabelPixelType label) const
1053 {
1054 MapConstIterator mapIt;
1055
1056 mapIt = m_LabelGeometryMapper.find(label);
1057 if ( mapIt == m_LabelGeometryMapper.end() )
1058 {
1059 // label does not exist, return a default value
1060 return nullptr;
1061 }
1062 else
1063 {
1064 return ( *mapIt ).second.m_OrientedIntensityImage;
1065 }
1066 }
1067
1068 template< typename TImage, typename TLabelImage >
1069 void
1070 LabelGeometryImageFilter< TImage, TLabelImage >
PrintSelf(std::ostream & os,Indent indent) const1071 ::PrintSelf(std::ostream & os, Indent indent) const
1072 {
1073 Superclass::PrintSelf(os, indent);
1074
1075 os << indent << "Number of labels: " << m_LabelGeometryMapper.size()
1076 << std::endl;
1077
1078 MapConstIterator mapIt;
1079 for ( mapIt = m_LabelGeometryMapper.begin(); mapIt != m_LabelGeometryMapper.end(); mapIt++ )
1080 {
1081 using LabelPrintType = typename NumericTraits< LabelPixelType >::PrintType;
1082 os << indent << "Label[" << (LabelPrintType)( ( *mapIt ).second.m_Label ) << "]: ";
1083 os << "\t Volume: " << ( *mapIt ).second.m_ZeroOrderMoment;
1084 os << "\t Integrated Intensity: " << ( *mapIt ).second.m_Sum;
1085 os << "\t Centroid: " << ( *mapIt ).second.m_Centroid;
1086 os << "\t Weighted Centroid: " << ( *mapIt ).second.m_WeightedCentroid;
1087 os << "\t Axes Length: " << ( *mapIt ).second.m_AxesLength;
1088 os << "\t Eccentricity: " << ( *mapIt ).second.m_Eccentricity;
1089 os << "\t Elongation: " << ( *mapIt ).second.m_Elongation;
1090 os << "\t Orientation: " << ( *mapIt ).second.m_Orientation;
1091 os << "\t Bounding box: " << ( *mapIt ).second.m_BoundingBox;
1092 os << "\t Bounding box volume: " << ( *mapIt ).second.m_BoundingBoxVolume;
1093 os << "\t Bounding box size: " << ( *mapIt ).second.m_BoundingBoxSize;
1094 // Oriented bounding box verticies
1095 os << "\t Oriented bounding box volume: " << ( *mapIt ).second.m_OrientedBoundingBoxVolume;
1096 os << "\t Oriented bounding box size: " << ( *mapIt ).second.m_OrientedBoundingBoxSize;
1097 // Rotation matrix
1098 os << std::endl;
1099 os << "\t Calculate oriented intensity regions: " << m_CalculateOrientedIntensityRegions;
1100 os << "\t Calculate pixel indices: " << m_CalculatePixelIndices;
1101 os << "\t Calculate oriented bounding box: " << m_CalculateOrientedBoundingBox;
1102 os << "\t Calculate oriented label regions: " << m_CalculateOrientedLabelRegions;
1103 os << "\n\n";
1104 }
1105 }
1106 } // end namespace itk
1107 #endif
1108