otbShapeAttributesLabelMapFilter.txx 41.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
/*=========================================================================

  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.

  Some parts of this code are derived from ITK. See ITKCopyright.txt
  for details.


Emmanuel Christophe's avatar
Emmanuel Christophe committed
16
17
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
18
19
20
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
21
22
#ifndef otbShapeAttributesLabelMapFilter_txx
#define otbShapeAttributesLabelMapFilter_txx
23
24
25

#include "otbShapeAttributesLabelMapFilter.h"
#include "itkProgressReporter.h"
26
27
#include "itkConstNeighborhoodIterator.h"
#include "itkConstShapedNeighborhoodIterator.h"
28
29
#include "itkLabelMapToLabelImageFilter.h"
#include "itkConstantBoundaryCondition.h"
30
31
#include "itkGeometryUtilities.h"
#include "itkConnectedComponentAlgorithm.h"
32
33
34
35
#include "vnl/algo/vnl_real_eigensystem.h"
#include "vnl/algo/vnl_symmetric_eigensystem.h"

#include "otbMacro.h"
36
#include <deque>
37
38
39

namespace otb {

40
41
42
namespace Functor {

template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
43
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
44
::ShapeAttributesLabelObjectFunctor() : m_ComputeFeretDiameter(false),
OTB Bot's avatar
STYLE    
OTB Bot committed
45
  m_ComputePerimeter(false),
46
  m_ComputeFlusser(true),
47
  m_ComputePolygon(true),
OTB Bot's avatar
STYLE    
OTB Bot committed
48
  m_ReducedAttributeSet(true),
49
  m_LabelImage(ITK_NULLPTR)
50
51
52
53
{}

/** The comparator (!=) */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
54
55
56
57
bool
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::operator != (const Self &self)
  {
58
59
60
61
62
63
64
65
66
67
68
  // Initialize response
  bool resp = true;

  // Check for differences
  resp = resp && (m_ComputeFeretDiameter != self.m_ComputeFeretDiameter);
  resp = resp && (m_ComputePerimeter != self.m_ComputePerimeter);
  resp = resp && (m_ReducedAttributeSet != self.m_ReducedAttributeSet);
  resp = resp && (m_LabelImage != self.m_LabelImage);

  // Return
  return resp;
OTB Bot's avatar
STYLE    
OTB Bot committed
69
  }
70

71
72
/** The comparator (==)*/
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
73
74
75
76
bool
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::operator == (const Self &self)
  {
77
78
  // Call the != implementation
  return !(this != self);
OTB Bot's avatar
STYLE    
OTB Bot committed
79
80
  }

81
82
/** Set the compute perimeter flag */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
83
84
void
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
85
86
87
88
::SetComputePerimeter(bool flag)
{
  m_ComputePerimeter = flag;
}
89

90
91
/** Get the compute perimeter flag */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
92
93
bool
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
94
95
96
97
98
::GetComputePerimeter() const
{
  return m_ComputePerimeter;
}

99
/** Set the compute polygon flag */
100
101
102
103
104
105
106
107
template <class TLabelObject, class TLabelImage>
void
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::SetComputePolygon(bool flag)
{
  m_ComputePolygon = flag;
}

108
/** Get the compute polygon flag */
109
110
111
112
113
114
115
116
template <class TLabelObject, class TLabelImage>
bool
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::GetComputePolygon() const
{
  return m_ComputePolygon;
}

117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
/** Set the compute Flusser flag */
template <class TLabelObject, class TLabelImage>
void
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::SetComputeFlusser(bool flag)
{
  m_ComputeFlusser = flag;
}

/** Get the compute Flusser flag */
template <class TLabelObject, class TLabelImage>
bool
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::GetComputeFlusser() const
{
  return m_ComputeFlusser;
}
134

135
136
/** Set the compute feret diameter flag */
template <class TLabelObject, class TLabelImage>
137
void
OTB Bot's avatar
STYLE    
OTB Bot committed
138
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
139
::SetComputeFeretDiameter(bool flag)
140
{
141
142
  m_ComputeFeretDiameter = flag;
}
OTB Bot's avatar
STYLE    
OTB Bot committed
143

144
145
/** Get the compute feret diameter flag */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
146
147
bool
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
148
149
150
151
::GetComputeFeretDiameter() const
{
  return m_ComputeFeretDiameter;
}
OTB Bot's avatar
STYLE    
OTB Bot committed
152

153
154
155
/** Set the compute reduced attributes set flag */
template <class TLabelObject, class TLabelImage>
void
OTB Bot's avatar
STYLE    
OTB Bot committed
156
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
157
158
159
160
::SetReducedAttributeSet(bool flag)
{
  m_ReducedAttributeSet = flag;
}
161

162
163
164
/** Get the compute reduced attributes set flag */
template <class TLabelObject, class TLabelImage>
bool
OTB Bot's avatar
STYLE    
OTB Bot committed
165
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
166
167
168
169
::GetReducedAttributeSet() const
{
  return m_ReducedAttributeSet;
}
170

OTB Bot's avatar
STYLE    
OTB Bot committed
171
/** Set the label image (used only to compute
172
173
174
 *  the Feret diameter */
template <class TLabelObject, class TLabelImage>
void
OTB Bot's avatar
STYLE    
OTB Bot committed
175
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
176
177
178
179
::SetLabelImage(const TLabelImage * image)
{
  m_LabelImage = image;
}
180

181
182
/** Get the label image */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
183
184
const TLabelImage *
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
185
186
187
::GetLabelImage() const
{
  return m_LabelImage;
188
189
190
}


191
/** This is the functor implementation
OTB Bot's avatar
STYLE    
OTB Bot committed
192
 *  Calling the functor on a label object
193
194
 *  will update its shape attributes */
template <class TLabelObject, class TLabelImage>
195
void
OTB Bot's avatar
STYLE    
OTB Bot committed
196
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
197
::operator() (LabelObjectType * lo)
198
{
OTB Bot's avatar
STYLE    
OTB Bot committed
199
  const typename LabelObjectType::LabelType& label = lo->GetLabel();
200
201
202
203
204

  // TODO: compute sizePerPixel, borderMin and borderMax in BeforeThreadedGenerateData() ?

  // compute the size per pixel, to be used later
  double sizePerPixel = 1;
205
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
206
    {
207
    sizePerPixel *= vcl_abs(m_LabelImage->GetSpacing()[i]);
208
    }
OTB Bot's avatar
STYLE    
OTB Bot committed
209
210

  typename std::vector<double> sizePerPixelPerDimension;
211
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
212
    {
OTB Bot's avatar
STYLE    
OTB Bot committed
213
    sizePerPixelPerDimension.push_back(sizePerPixel / vcl_abs(m_LabelImage->GetSpacing()[i]));
214
    }
OTB Bot's avatar
STYLE    
OTB Bot committed
215

216
  // compute the max the index on the border of the image
217
218
  typename LabelObjectType::IndexType borderMin = m_LabelImage->GetLargestPossibleRegion().GetIndex();
  typename LabelObjectType::IndexType borderMax = borderMin;
219
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
220
    {
221
    borderMax[i] += borderMin[i] + m_LabelImage->GetLargestPossibleRegion().GetSize()[i] - 1;
222
223
224
    }

  // init the vars
OTB Bot's avatar
STYLE    
OTB Bot committed
225
226
227
  unsigned long                                                 size = 0;
  itk::ContinuousIndex<double, LabelObjectType::ImageDimension> centroid;
  centroid.Fill(0);
228
  typename LabelObjectType::IndexType mins;
OTB Bot's avatar
STYLE    
OTB Bot committed
229
  mins.Fill(itk::NumericTraits<long>::max());
230
  typename LabelObjectType::IndexType maxs;
OTB Bot's avatar
STYLE    
OTB Bot committed
231
232
233
234
235
  maxs.Fill(itk::NumericTraits<long>::NonpositiveMin());
  unsigned long                                                                         sizeOnBorder = 0;
  double                                                                                physicalSizeOnBorder = 0;
  itk::Matrix<double, LabelObjectType::ImageDimension, LabelObjectType::ImageDimension> centralMoments;
  centralMoments.Fill(0);
236

237
238
239
240
  ConstLineIteratorType lit = ConstLineIteratorType(lo);
  lit.GoToBegin();
  //lit = IteratorType( ref );
  //typename LabelObjectType::LineContainerType&                lineContainer = lo->GetLineContainer();
241
242

  // iterate over all the lines
243
244
  while ( !lit.IsAtEnd() )
    //for (lit = lineContainer.begin(); lit != lineContainer.end(); lit++)
245
    {
246
247
    const typename LabelObjectType::IndexType& idx = lit.GetLine().GetIndex();
    unsigned long                              length = lit.GetLine().GetLength();
248
249
250
251
252
253

    // update the size
    size += length;

    // update the centroid - and report the progress
    // first, update the axes which are not 0
254
    for (DimensionType i = 1; i < LabelObjectType::ImageDimension; ++i)
255
256
257
258
      {
      centroid[i] += length * idx[i];
      }
    // then, update the axis 0
OTB Bot's avatar
STYLE    
OTB Bot committed
259
    centroid[0] += idx[0] * length + (length * (length - 1)) / 2.0;
260
261

    // update the mins and maxs
262
    for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
263
      {
OTB Bot's avatar
STYLE    
OTB Bot committed
264
      if (idx[i] < mins[i])
265
266
267
        {
        mins[i] = idx[i];
        }
OTB Bot's avatar
STYLE    
OTB Bot committed
268
      if (idx[i] > maxs[i])
269
270
271
272
273
        {
        maxs[i] = idx[i];
        }
      }
    // must fix the max for the axis 0
OTB Bot's avatar
STYLE    
OTB Bot committed
274
    if (idx[0] + (long) length > maxs[0])
275
276
277
278
279
280
      {
      maxs[0] = idx[0] + length - 1;
      }

    // object is on a border ?
    bool isOnBorder = false;
281
    for (DimensionType i = 1; i < LabelObjectType::ImageDimension; ++i)
282
      {
OTB Bot's avatar
STYLE    
OTB Bot committed
283
      if (idx[i] == borderMin[i] || idx[i] == borderMax[i])
284
285
286
287
288
        {
        isOnBorder = true;
        break;
        }
      }
OTB Bot's avatar
STYLE    
OTB Bot committed
289
    if (isOnBorder)
290
291
292
293
294
295
296
297
298
      {
      // the line touch a border on a dimension other than 0, so
      // all the line touch a border
      sizeOnBorder += length;
      }
    else
      {
      // we must check for the dimension 0
      bool isOnBorder0 = false;
OTB Bot's avatar
STYLE    
OTB Bot committed
299
      if (idx[0] == borderMin[0])
300
301
302
303
304
        {
        // one more pixel on the border
        sizeOnBorder++;
        isOnBorder0 = true;
        }
OTB Bot's avatar
STYLE    
OTB Bot committed
305
      if (!isOnBorder0 || length > 1)
306
307
        {
        // we can check for the end of the line
OTB Bot's avatar
STYLE    
OTB Bot committed
308
        if (idx[0] + (long) length - 1 == borderMax[0])
309
310
311
312
313
314
          {
          // one more pixel on the border
          sizeOnBorder++;
          }
        }
      }
OTB Bot's avatar
STYLE    
OTB Bot committed
315

316
317
    // physical size on border
    // first, the dimension 0
OTB Bot's avatar
STYLE    
OTB Bot committed
318
    if (idx[0] == borderMin[0])
319
320
321
322
      {
      // the begining of the line
      physicalSizeOnBorder += sizePerPixelPerDimension[0];
      }
OTB Bot's avatar
STYLE    
OTB Bot committed
323
    if (idx[0] + (long) length - 1 == borderMax[0])
324
325
326
327
328
      {
      // and the end of the line
      physicalSizeOnBorder += sizePerPixelPerDimension[0];
      }
    // then the other dimensions
329
    for (DimensionType i = 1; i < LabelObjectType::ImageDimension; ++i)
330
      {
OTB Bot's avatar
STYLE    
OTB Bot committed
331
      if (idx[i] == borderMin[i])
332
333
334
335
        {
        // one border
        physicalSizeOnBorder += sizePerPixelPerDimension[i] * length;
        }
OTB Bot's avatar
STYLE    
OTB Bot committed
336
      if (idx[i] == borderMax[i])
337
338
339
340
341
        {
        // and the other
        physicalSizeOnBorder += sizePerPixelPerDimension[i] * length;
        }
      }
342
343

    /*
344
    // moments computation
345
    // ****************************************************************
346
347
    // that commented code is the basic implementation. The next piece of code
    // give the same result in a much efficient way, by using expanded formulae
348
349
    // allowed by the binary case instead of loops.
    // ****************************************************************
350
351
352
353
354
355
         long endIdx0 = idx[0] + length;
         for( typename LabelObjectType::IndexType iidx = idx; iidx[0]<endIdx0; iidx[0]++)
           {
           typename TLabelImage::PointType pP;
           m_LabelImage->TransformIndexToPhysicalPoint(iidx, pP);

356
           for(unsigned int i=0; i<LabelObjectType::ImageDimension; ++i)
357
             {
358
             for(unsigned int j=0; j<LabelObjectType::ImageDimension; ++j)
359
360
361
362
363
364
365
               {
               centralMoments[i][j] += pP[i] * pP[j];
               }
             }
           }
    */

366
    // get the physical position and the spacing - they are used several times later
367
    typename TLabelImage::PointType physicalPosition;
OTB Bot's avatar
STYLE    
OTB Bot committed
368
369
    m_LabelImage->TransformIndexToPhysicalPoint(idx, physicalPosition);
    const typename TLabelImage::SpacingType& spacing = m_LabelImage->GetSpacing();
370
    // the sum of x positions, also reused several times
OTB Bot's avatar
STYLE    
OTB Bot committed
371
    double sumX = length * (physicalPosition[0] + (spacing[0] * (length - 1)) / 2.0);
372
373
    // the real job - the sum of square of x positions
    // that's the central moments for dims 0, 0
OTB Bot's avatar
STYLE    
OTB Bot committed
374
375
376
    centralMoments[0][0] += length * (physicalPosition[0] * physicalPosition[0]
                                      + spacing[0] *
                                      (length - 1) * ((spacing[0] * (2 * length - 1)) / 6.0 + physicalPosition[0]));
377
    // the other ones
378
    for (DimensionType i = 1; i < LabelObjectType::ImageDimension; ++i)
379
380
381
382
      {
      // do this one here to avoid the double assigment in the following loop
      // when i == j
      centralMoments[i][i] += length * physicalPosition[i] * physicalPosition[i];
383
      // central moments are symmetrics, so avoid to compute them 2 times
384
      for (DimensionType j = i + 1; j < LabelObjectType::ImageDimension; ++j)
385
386
387
388
389
390
391
392
393
394
395
396
        {
        // note that we won't use that code if the image dimension is less than 3
        // --> the tests should be in 3D at least
        double cm = length * physicalPosition[i] * physicalPosition[j];
        centralMoments[i][j] += cm;
        centralMoments[j][i] += cm;
        }
      // the last moments: the ones for the dimension 0
      double cm = sumX * physicalPosition[i];
      centralMoments[i][0] += cm;
      centralMoments[0][i] += cm;
      }
397
    ++lit;
398
399
400
    }

  // final computation
401
  typename TLabelImage::SizeType regionSize;
OTB Bot's avatar
STYLE    
OTB Bot committed
402
403
  double                         minSize = itk::NumericTraits<double>::max();
  double                         maxSize = itk::NumericTraits<double>::NonpositiveMin();
404
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
405
406
407
    {
    centroid[i] /= size;
    regionSize[i] = maxs[i] - mins[i] + 1;
408
    double s = regionSize[i] * vcl_abs(m_LabelImage->GetSpacing()[i]);
OTB Bot's avatar
STYLE    
OTB Bot committed
409
410
    minSize = std::min(s, minSize);
    maxSize = std::max(s, maxSize);
411
    for (DimensionType j = 0; j < LabelObjectType::ImageDimension; ++j)
412
413
414
415
      {
      centralMoments[i][j] /= size;
      }
    }
OTB Bot's avatar
STYLE    
OTB Bot committed
416
  typename TLabelImage::RegionType region(mins, regionSize);
417
  typename TLabelImage::PointType physicalCentroid;
OTB Bot's avatar
STYLE    
OTB Bot committed
418
  m_LabelImage->TransformContinuousIndexToPhysicalPoint(centroid, physicalCentroid);
419
420

  // Center the second order moments
421
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
422
    {
423
    for (DimensionType j = 0; j < LabelObjectType::ImageDimension; ++j)
424
425
426
427
428
429
      {
      centralMoments[i][j] -= physicalCentroid[i] * physicalCentroid[j];
      }
    }

  // Compute principal moments and axes
OTB Bot's avatar
STYLE    
OTB Bot committed
430
431
  itk::Vector<double, LabelObjectType::ImageDimension> principalMoments;
  vnl_symmetric_eigensystem<double> eigen(centralMoments.GetVnlMatrix());
432
  vnl_diag_matrix<double> pm = eigen.D;
433
  for (unsigned int i = 0; i < LabelObjectType::ImageDimension; ++i)
434
    {
OTB Bot's avatar
STYLE    
OTB Bot committed
435
//    principalMoments[i] = 4 * vcl_sqrt( pm(i, i) );
OTB Bot's avatar
STYLE    
OTB Bot committed
436
    principalMoments[i] = pm(i, i);
437
    }
OTB Bot's avatar
STYLE    
OTB Bot committed
438
439
  itk::Matrix<double, LabelObjectType::ImageDimension, LabelObjectType::ImageDimension>
  principalAxes = eigen.V.transpose();
440
441
442

  // Add a final reflection if needed for a proper rotation,
  // by multiplying the last row by the determinant
OTB Bot's avatar
STYLE    
OTB Bot committed
443
444
445
  vnl_real_eigensystem eigenrot(principalAxes.GetVnlMatrix());
  vnl_diag_matrix<vcl_complex<double> > eigenval = eigenrot.D;
  vcl_complex<double> det(1.0, 0.0);
446

447
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
448
    {
OTB Bot's avatar
STYLE    
OTB Bot committed
449
    det *= eigenval(i, i);
450
451
    }

452
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
453
    {
OTB Bot's avatar
STYLE    
OTB Bot committed
454
    principalAxes[LabelObjectType::ImageDimension - 1][i] *= std::real(det);
455
456
457
    }

  double elongation = 0;
OTB Bot's avatar
STYLE    
OTB Bot committed
458
  if (principalMoments[LabelObjectType::ImageDimension - 2] != 0)
459
    {
OTB Bot's avatar
STYLE    
OTB Bot committed
460
461
462
    elongation =
      vcl_sqrt(principalMoments[LabelObjectType::ImageDimension -
                                1] / principalMoments[LabelObjectType::ImageDimension - 2]);
463
464
465
    }

  double physicalSize = size * sizePerPixel;
OTB Bot's avatar
STYLE    
OTB Bot committed
466
467
  double equivalentRadius = hyperSphereRadiusFromVolume(physicalSize);
  double equivalentPerimeter = hyperSpherePerimeter(equivalentRadius);
468
469

  // compute equilalent ellipsoid radius
OTB Bot's avatar
STYLE    
OTB Bot committed
470
471
  itk::Vector<double, LabelObjectType::ImageDimension> ellipsoidSize;
  double                                               edet = 1.0;
472
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
473
474
475
    {
    edet *= principalMoments[i];
    }
OTB Bot's avatar
STYLE    
OTB Bot committed
476
  edet = vcl_pow(edet, 1.0 / LabelObjectType::ImageDimension);
477
  for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
478
    {
OTB Bot's avatar
STYLE    
OTB Bot committed
479
    ellipsoidSize[i] = 2.0 * equivalentRadius * vcl_sqrt(principalMoments[i] / edet);
480
481
    }

482

483
484
485
486
487
488
489
490
  if (m_ComputeFlusser)
  {
    // Flusser moments (only make sense when ImageDimension == 2)
    if (LabelObjectType::ImageDimension == 2)
      {
      // complex centered moments
      std::complex<double> c11, c20, c12, c21, c22, c30, c31, c40;
      c11 = c20 = c12 = c21 = c22 = c30 = c31 = c40 = std::complex<double>(0, 0);
491
492
493
494
495
      // Update the line iterator to the beginning
      lit.GoToBegin();
      //for (lit = lineContainer.begin(); lit != lineContainer.end();
      //lit++)
      while ( ! lit.IsAtEnd() )
496
        {
497
498
        const typename LabelObjectType::IndexType& idx = lit.GetLine().GetIndex();
        unsigned long length = lit.GetLine().GetLength();
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
        //
        long endIdx0 = idx[0] + length;
        for (typename LabelObjectType::IndexType iidx = idx; iidx[0] < endIdx0; iidx[0]++)
          {
          typename TLabelImage::PointType cPP;
          m_LabelImage->TransformIndexToPhysicalPoint(iidx, cPP);
          cPP -= physicalCentroid.GetVectorFromOrigin();
          std::complex<double> xpiy(cPP[0], cPP[1]); // x + i y
          std::complex<double> xmiy(cPP[0], -cPP[1]); // x - i y

          c11 += xpiy * xmiy;
          c20 += xpiy * xpiy;
          c12 += xpiy * xmiy * xmiy;
          c21 += xpiy * xpiy * xmiy;
          c30 += xpiy * xpiy * xpiy;
          c22 += xpiy * xpiy * xmiy * xmiy;
          c31 += xpiy * xpiy * xpiy * xmiy;
          c40 += xpiy * xpiy * xpiy * xpiy;
          }
OTB Bot's avatar
STYLE    
OTB Bot committed
518
       ++lit;
519
520
        }

521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
      // normalize
      c11 /= physicalSize * physicalSize;
      c20 /= physicalSize * physicalSize;
      c12 /= vcl_pow(physicalSize, 5.0 / 2);
      c21 /= vcl_pow(physicalSize, 5.0 / 2);
      c30 /= vcl_pow(physicalSize, 5.0 / 2);
      c22 /= vcl_pow(physicalSize, 3);
      c31 /= vcl_pow(physicalSize, 3);
      c40 /= vcl_pow(physicalSize, 3);

      lo->SetAttribute("SHAPE::Flusser01", c11.real());
      lo->SetAttribute("SHAPE::Flusser02", (c21 * c12).real());
      lo->SetAttribute("SHAPE::Flusser03", (c20 * vcl_pow(c12, 2)).real());
      lo->SetAttribute("SHAPE::Flusser04", (c20 * vcl_pow(c12, 2)).imag());
      lo->SetAttribute("SHAPE::Flusser05", (c30 * vcl_pow(c12, 3)).real());
      lo->SetAttribute("SHAPE::Flusser06", (c30 * vcl_pow(c12, 3)).imag());
      lo->SetAttribute("SHAPE::Flusser07", c22.real());
      lo->SetAttribute("SHAPE::Flusser08", (c31 * vcl_pow(c12, 2)).real());
      lo->SetAttribute("SHAPE::Flusser09", (c31 * vcl_pow(c12, 2)).imag());
      lo->SetAttribute("SHAPE::Flusser10", (c40 * vcl_pow(c12, 4)).real());
      lo->SetAttribute("SHAPE::Flusser11", (c40 * vcl_pow(c12, 4)).imag());
      }
543
    }
544

545
  // Set the attributes
546
547
548
549
550
551
552
553
554
555
  if (m_ComputePolygon)
    {
    PolygonFunctorType polygonFunctor;
    SimplifyPolygonFunctorType simplifyFunctor;
    polygonFunctor.SetStartIndex(m_LabelImage->GetLargestPossibleRegion().GetIndex());
    polygonFunctor.SetOrigin(m_LabelImage->GetOrigin());
    polygonFunctor.SetSpacing(m_LabelImage->GetSpacing());
    typename PolygonType::Pointer polygon = simplifyFunctor(polygonFunctor(lo));
    lo->SetPolygon(polygon);
    }
556
557

  // Physical size
OTB Bot's avatar
STYLE    
OTB Bot committed
558
559
  lo->SetAttribute("SHAPE::PhysicalSize", physicalSize);

560
  // Elongation
OTB Bot's avatar
STYLE    
OTB Bot committed
561
  lo->SetAttribute("SHAPE::Elongation", elongation);
562

OTB Bot's avatar
STYLE    
OTB Bot committed
563
  if (m_ComputeFeretDiameter)
564
565
    {
    // init the vars
566
    unsigned long ssize = 0;
OTB Bot's avatar
STYLE    
OTB Bot committed
567
    typedef typename std::deque<typename LabelObjectType::IndexType> IndexListType;
568
    IndexListType idxList;
OTB Bot's avatar
STYLE    
OTB Bot committed
569

570
    // Line iterator
571
572
    ConstLineIteratorType llit = ConstLineIteratorType(lo);
    llit.GoToBegin();
573

OTB Bot's avatar
STYLE    
OTB Bot committed
574
    typedef typename itk::ConstNeighborhoodIterator<LabelImageType> NeighborIteratorType;
575
    typename TLabelImage::SizeType neighborHoodRadius;
OTB Bot's avatar
STYLE    
OTB Bot committed
576
577
    neighborHoodRadius.Fill(1);
    NeighborIteratorType it(neighborHoodRadius, m_LabelImage, m_LabelImage->GetBufferedRegion());
578
579
    itk::ConstantBoundaryCondition<LabelImageType> lcbc;
    // use label + 1 to have a label different of the current label on the border
OTB Bot's avatar
STYLE    
OTB Bot committed
580
581
    lcbc.SetConstant(label + 1);
    it.OverrideBoundaryCondition(&lcbc);
582
583
584
    it.GoToBegin();

    // iterate over all the lines
585
    while ( !llit.IsAtEnd() )
586
      {
587
588
      const typename LabelObjectType::IndexType& firstIdx = llit.GetLine().GetIndex();
      unsigned long                              length = llit.GetLine().GetLength();
589
590

      long endIdx0 = firstIdx[0] + length;
OTB Bot's avatar
STYLE    
OTB Bot committed
591
      for (typename LabelObjectType::IndexType idx = firstIdx; idx[0] < endIdx0; idx[0]++)
592
593
594
595
596
597
        {

        // move the iterator to the new location
        it += idx - it.GetIndex();

        // push the pixel in the list if it is on the border of the object
598
        for (unsigned i = 0; i < it.Size(); ++i)
599
          {
OTB Bot's avatar
STYLE    
OTB Bot committed
600
          if (it.GetPixel(i) != label)
601
            {
OTB Bot's avatar
STYLE    
OTB Bot committed
602
            idxList.push_back(idx);
603
            ssize++;
604
605
606
607
            break;
            }
          }
        }
608
      ++llit;
609
610
611
612
      }

    // we can now search the feret diameter
    double feretDiameter = 0;
OTB Bot's avatar
STYLE    
OTB Bot committed
613
614
615
    for (typename IndexListType::const_iterator iIt1 = idxList.begin();
         iIt1 != idxList.end();
         iIt1++)
616
617
      {
      typename IndexListType::const_iterator iIt2 = iIt1;
OTB Bot's avatar
STYLE    
OTB Bot committed
618
      for (iIt2++; iIt2 != idxList.end(); iIt2++)
619
620
621
        {
        // Compute the length between the 2 indexes
        double length = 0;
622
        for (DimensionType i = 0; i < LabelObjectType::ImageDimension; ++i)
623
          {
OTB Bot's avatar
STYLE    
OTB Bot committed
624
          length += vcl_pow((iIt1->operator[] (i) - iIt2->operator[] (i)) * m_LabelImage->GetSpacing()[i], 2);
625
          }
OTB Bot's avatar
STYLE    
OTB Bot committed
626
        if (feretDiameter < length)
627
628
629
630
631
632
          {
          feretDiameter = length;
          }
        }
      }
    // final computation
OTB Bot's avatar
STYLE    
OTB Bot committed
633
    feretDiameter = vcl_sqrt(feretDiameter);
634
635

    // finally put the values in the label object
OTB Bot's avatar
STYLE    
OTB Bot committed
636
    lo->SetAttribute("SHAPE::FeretDiameter", feretDiameter);
637
638
    }

639
  // be sure that the calculator has the perimeter estimation for that label.
640
  // The calculator may not have the label if the object is only on a border.
641
  // It will occur for sure when processing a 2D image with a 3D filter.
642
  if (m_ComputePerimeter)
643
    {
644
645
    double perimeter = this->ComputePerimeter(lo,region);
    //double perimeter = lo->ComputePerimeter();
OTB Bot's avatar
STYLE    
OTB Bot committed
646
    lo->SetAttribute("SHAPE::Perimeter", perimeter);
647
    lo->SetAttribute("SHAPE::Roundness", equivalentPerimeter / perimeter);
648
649
650
651
    }

  // Complete feature set

OTB Bot's avatar
STYLE    
OTB Bot committed
652
  if (!m_ReducedAttributeSet)
653
    {
OTB Bot's avatar
STYLE    
OTB Bot committed
654
    lo->SetAttribute("SHAPE::Size", size);
Emmanuel Christophe's avatar
Emmanuel Christophe committed
655
    std::ostringstream oss;
OTB Bot's avatar
STYLE    
OTB Bot committed
656
    for (unsigned int dim = 0; dim < LabelObjectType::ImageDimension; ++dim)
657
658
      {
      oss.str("");
OTB Bot's avatar
STYLE    
OTB Bot committed
659
      oss << "SHAPE::RegionIndex" << dim;
660
      lo->SetAttribute(oss.str().c_str(), region.GetIndex()[dim]);
661
      oss.str("");
OTB Bot's avatar
STYLE    
OTB Bot committed
662
      oss << "SHAPE::RegionSize" << dim;
663
      lo->SetAttribute(oss.str().c_str(), region.GetSize()[dim]);
664
      oss.str("");
OTB Bot's avatar
STYLE    
OTB Bot committed
665
      oss << "SHAPE::PhysicalCentroid" << dim;
666
      lo->SetAttribute(oss.str().c_str(), physicalCentroid[dim]);
667
      oss.str("");
OTB Bot's avatar
STYLE    
OTB Bot committed
668
669
      oss << "SHAPE::EquivalentEllipsoidRadius" << dim;
      lo->SetAttribute(oss.str().c_str(), ellipsoidSize[dim]);
670
      oss.str("");
OTB Bot's avatar
STYLE    
OTB Bot committed
671
672
673
674
675
676
677
678
679
      oss << "SHAPE::PrincipalMoments" << dim;
      lo->SetAttribute(oss.str().c_str(), principalMoments[dim]);

      for (unsigned int dim2 = 0; dim2 < LabelObjectType::ImageDimension; ++dim2)
        {
        oss.str("");
        oss << "SHAPE::PrincipalAxis" << dim << dim2;
        lo->SetAttribute(oss.str().c_str(), principalAxes(dim, dim2));
        }
680
      }
OTB Bot's avatar
STYLE    
OTB Bot committed
681
682
683
684
685
686

    lo->SetAttribute("SHAPE::RegionElongation", maxSize / minSize);
    lo->SetAttribute("SHAPE::RegionRatio", size / (double) region.GetNumberOfPixels());
    lo->SetAttribute("SHAPE::SizeOnBorder", sizeOnBorder);
    lo->SetAttribute("SHAPE::PhysicalSizeOnBorder", physicalSizeOnBorder);
    lo->SetAttribute("SHAPE::EquivalentPerimeter", equivalentPerimeter);
687
    lo->SetAttribute("SHAPE::EquivalentRadius",    equivalentRadius);
688
689
690
    }
}

691
692
693
694
695
696
697
698
699
700
701
702
template <class TLabelObject, class TLabelImage>
double
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::ComputePerimeter(LabelObjectType *labelObject, const RegionType & region)
{
  // store the lines in a N-1D image of vectors
  typedef std::deque< typename LabelObjectType::LineType > VectorLineType;
  typedef itk::Image< VectorLineType, ImageDimension - 1 > LineImageType;
  typename LineImageType::Pointer lineImage = LineImageType::New();
  typename LineImageType::IndexType lIdx;
  typename LineImageType::SizeType lSize;
  RegionType boundingBox = region;
703
  for( DimensionType i=0; i<ImageDimension-1; i++ )
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
    {
    lIdx[i] = boundingBox.GetIndex()[i+1];
    lSize[i] = boundingBox.GetSize()[i+1];
    }
  typename LineImageType::RegionType lRegion;
  lRegion.SetIndex( lIdx );
  lRegion.SetSize( lSize );
  // enlarge the region a bit to avoid boundary problems
  typename LineImageType::RegionType elRegion(lRegion);
  lSize.Fill(1);
  elRegion.PadByRadius(lSize);
  // std::cout << boundingBox << "  " << lRegion << "  " << elRegion << std::endl;
  // now initialize the image
  lineImage->SetRegions( elRegion );
  lineImage->Allocate();
  lineImage->FillBuffer( VectorLineType() );

  // std::cout << "lineContainer.size(): " << lineContainer.size() << std::endl;

  // Iterate over all the lines and fill the image of lines
  typename LabelObjectType::ConstLineIterator lit( labelObject );
  while( ! lit.IsAtEnd() )
    {
    const typename TLabelObject::IndexType & idx = lit.GetLine().GetIndex();
728
    for( DimensionType i=0; i<ImageDimension-1; i++ )
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
      {
      lIdx[i] = idx[i+1];
      }
    lineImage->GetPixel( lIdx ).push_back( lit.GetLine() );
    ++lit;
    }

  // a data structure to store the number of intercepts on each direction
  typedef typename std::map<OffsetType, itk::SizeValueType, typename OffsetType::LexicographicCompare> MapInterceptType;
  MapInterceptType intercepts;
  // int nbOfDirections = (int)vcl_pow( 2.0, (int)ImageDimension ) - 1;
  // intecepts.resize(nbOfDirections + 1);  // code begins at position 1

  // now iterate over the vectors of lines
  typedef itk::ConstShapedNeighborhoodIterator< LineImageType > LineImageIteratorType;
  LineImageIteratorType lIt( lSize, lineImage, lRegion ); // the original, non padded region
  setConnectivity( &lIt, true );
  for( lIt.GoToBegin(); !lIt.IsAtEnd(); ++lIt )
    {
    const VectorLineType & ls = lIt.GetCenterPixel();

    // there are two intercepts on the 0 axis for each line
    OffsetType no;
    no.Fill(0);
    no[0] = 1;
    // std::cout << no << "-> " << 2 * ls.size() << std::endl;
    intercepts[no] += 2 * ls.size();

    // and look at the neighbors
    typename LineImageIteratorType::ConstIterator ci;
    for (ci = lIt.Begin(); ci != lIt.End(); ci++)
      {
          // std::cout << "-------------" << std::endl;
      // the vector of lines in the neighbor
      const VectorLineType & ns = ci.Get();
      // prepare the offset to be stored in the intercepts map
      typename LineImageType::OffsetType lno = ci.GetNeighborhoodOffset();
      no[0] = 0;
767
      for( DimensionType i=0; i<ImageDimension-1; i++ )
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
        {
        no[i+1] = vnl_math_abs(lno[i]);
        }
      OffsetType dno = no; // offset for the diagonal
      dno[0] = 1;

      // now process the two lines to search the pixels on the contour of the object
      if( ls.empty() )
        {
        // std::cout << "ls.empty()" << std::endl;
        // nothing to do
        }
      if( ns.empty() )
        {
        // no line in the neighbors - all the lines in ls are on the contour
        for( typename VectorLineType::const_iterator li = ls.begin(); li != ls.end(); ++li )
          {
          // std::cout << "ns.empty()" << std::endl;
          const typename LabelObjectType::LineType & l = *li;
          // add as much intercepts as the line size
          intercepts[no] += l.GetLength();
          // and 2 times as much diagonal intercepts as the line size
          intercepts[dno] += l.GetLength() * 2;
          }
        }
      else
        {
        // std::cout << "else" << std::endl;
        // TODO - fix the code when the line starts at  NumericTraits<IndexValueType>::NonpositiveMin()
        // or end at  NumericTraits<IndexValueType>::max()
        typename VectorLineType::const_iterator li = ls.begin();
        typename VectorLineType::const_iterator ni = ns.begin();

        itk::IndexValueType lZero = 0;
        itk::IndexValueType lMin = 0;
        itk::IndexValueType lMax = 0;

        itk::IndexValueType nMin = itk::NumericTraits<itk::IndexValueType>::NonpositiveMin() + 1;
        itk::IndexValueType nMax = ni->GetIndex()[0] - 1;

        while( li!=ls.end() )
          {
          // update the current line min and max. Neighbor line data is already up to date.
          lMin = li->GetIndex()[0];
          lMax = lMin + li->GetLength() - 1;

          // add as much intercepts as intersections of the 2 lines
          intercepts[no] += vnl_math_max( lZero, vnl_math_min(lMax, nMax) - vnl_math_max(lMin, nMin) + 1 );
          // std::cout << "============" << std::endl;
          // std::cout << "  lMin:" << lMin << " lMax:" << lMax << " nMin:" << nMin << " nMax:" << nMax;
          // std::cout << " count: " << vnl_math_max( 0l, vnl_math_min(lMax, nMax) - vnl_math_max(lMin, nMin) + 1 ) << std::endl;
          // std::cout << "  " << no << ": " << intercepts[no] << std::endl;
          // std::cout << vnl_math_max( lZero, vnl_math_min(lMax, nMax+1) - vnl_math_max(lMin, nMin+1) + 1 ) << std::endl;
          // std::cout << vnl_math_max( lZero, vnl_math_min(lMax, nMax-1) - vnl_math_max(lMin, nMin-1) + 1 ) << std::endl;
          // left diagonal intercepts
          intercepts[dno] += vnl_math_max( lZero, vnl_math_min(lMax, nMax+1) - vnl_math_max(lMin, nMin+1) + 1 );
          // right diagonal intercepts
          intercepts[dno] += vnl_math_max( lZero, vnl_math_min(lMax, nMax-1) - vnl_math_max(lMin, nMin-1) + 1 );

          // go to the next line or the next neighbor depending on where we are
          if(nMax <= lMax )
            {
            // go to next neighbor
            nMin = ni->GetIndex()[0] + ni->GetLength();
            ni++;

            if( ni != ns.end() )
              {
              nMax = ni->GetIndex()[0] - 1;
              }
            else
              {
              nMax = itk::NumericTraits<itk::IndexValueType>::max() - 1;
              }
            }
          else
            {
            // go to next line
            li++;
            }
          }

        }
      }
    }

  // compute the perimeter based on the intercept counts
  double perimeter = PerimeterFromInterceptCount( intercepts, m_LabelImage->GetSpacing() );
  return perimeter;
}

template <class TLabelObject, class TLabelImage>
template<class TMapIntercept, class TSpacing>
double
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::PerimeterFromInterceptCount( TMapIntercept & intercepts, const TSpacing & spacing )
{
  // std::cout << "PerimeterFromInterceptCount<>" << std::endl;
  double perimeter = 0.0;
  double pixelSize = 1.0;
  int dim = TSpacing::GetVectorDimension();
  for ( int i = 0; i < dim; i++ )
    {
    pixelSize *= spacing[i];
    }

  for( int i=0; i<dim; i++ )
    {
    OffsetType no;
    no.Fill(0);
    no[i] = 1;
    // std::cout << no << ": " << intercepts[no] << std::endl;
    perimeter += pixelSize / spacing[i] * intercepts[no]/2.0;
    }

  // Crofton's constant
  perimeter *= itk::GeometryUtilities::HyperSphereVolume( dim, 1.0 )
                 / itk::GeometryUtilities::HyperSphereVolume( dim - 1, 1.0 );
  return perimeter;
}

#if ! defined(ITK_DO_NOT_USE_PERIMETER_SPECIALIZATION)
template <class TLabelObject, class TLabelImage>
double
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::PerimeterFromInterceptCount( MapIntercept2Type & intercepts, const Spacing2Type spacing )
{
  // std::cout << "PerimeterFromInterceptCount2" << std::endl;
  double dx = spacing[0];
  double dy = spacing[1];

  Offset2Type nx =  {{1, 0}};
  Offset2Type ny =  {{0, 1}};
  Offset2Type nxy = {{1, 1}};

  // std::cout << "nx: " << intercepts[nx] << std::endl;
  // std::cout << "ny: " << intercepts[ny] << std::endl;
  // std::cout << "nxy: " << intercepts[nxy] << std::endl;

  double perimeter = 0.0;
  perimeter += dy * intercepts[nx]/2.0;
  perimeter += dx * intercepts[ny]/2.0;
  perimeter += dx*dy / spacing.GetNorm() * intercepts[nxy]/2.0;
  perimeter *= itk::Math::pi / 4.0;
  return perimeter;
};

template <class TLabelObject, class TLabelImage>
double
ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::PerimeterFromInterceptCount( MapIntercept3Type & intercepts, const Spacing3Type spacing )
{
  // std::cout << "PerimeterFromInterceptCount3" << std::endl;
  double dx = spacing[0];
  double dy = spacing[1];
  double dz = spacing[2];
  double dxy = vcl_sqrt( spacing[0]*spacing[0] + spacing[1]*spacing[1] );
  double dxz = vcl_sqrt( spacing[0]*spacing[0] + spacing[2]*spacing[2] );
  double dyz = vcl_sqrt( spacing[1]*spacing[1] + spacing[2]*spacing[2] );
  double dxyz = vcl_sqrt( spacing[0]*spacing[0] + spacing[1]*spacing[1] + spacing[2]*spacing[2] );
  double vol = spacing[0]*spacing[1]*spacing[2];

  // 'magical numbers', corresponding to area of voronoi partition on the
  // unit sphere, when germs are the 26 directions on the unit cube
  // Sum of (c1+c2+c3 + c4*2+c5*2+c6*2 + c7*4) equals 1.
  double c1 = 0.04577789120476 * 2;  // Ox
  double c2 = 0.04577789120476 * 2;  // Oy
  double c3 = 0.04577789120476 * 2;  // Oz
  double c4 = 0.03698062787608 * 2;  // Oxy
  double c5 = 0.03698062787608 * 2;  // Oxz
  double c6 = 0.03698062787608 * 2;  // Oyz
  double c7 = 0.03519563978232 * 2;  // Oxyz
  // TODO - recompute those values if the spacing is non isotrope

  Offset3Type nx =   {{1, 0, 0}};
  Offset3Type ny =   {{0, 1, 0}};
  Offset3Type nz =   {{0, 0, 1}};
  Offset3Type nxy =  {{1, 1, 0}};
  Offset3Type nxz =  {{1, 0, 1}};
  Offset3Type nyz =  {{0, 1, 1}};
  Offset3Type nxyz = {{1, 1, 1}};

  // std::cout << "nx: " << intercepts[nx] << std::endl;
  // std::cout << "ny: " << intercepts[ny] << std::endl;
  // std::cout << "nz: " << intercepts[nz] << std::endl;
  // std::cout << "nxy: " << intercepts[nxy] << std::endl;
  // std::cout << "nxz: " << intercepts[nxz] << std::endl;
  // std::cout << "nyz: " << intercepts[nyz] << std::endl;
  // std::cout << "nxyz: " << intercepts[nxyz] << std::endl;

  double perimeter = 0.0;
  perimeter += vol/dx * intercepts[nx]/2.0 * c1;
  perimeter += vol/dy * intercepts[ny]/2.0 * c2;
  perimeter += vol/dz * intercepts[nz]/2.0 * c3;
  perimeter += vol/dxy * intercepts[nxy]/2.0 * c4;
  perimeter += vol/dxz * intercepts[nxz]/2.0 * c5;
  perimeter += vol/dyz * intercepts[nyz]/2.0 * c6;
  perimeter += vol/dxyz * intercepts[nxyz]/2.0 * c7;
  perimeter *= 4;
  return perimeter;
};
#endif

971
972
/** Convenience internal method */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
973
974
long ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::factorial(long n)
975
{
OTB Bot's avatar
STYLE    
OTB Bot committed
976
  if (n < 1)
977
978
979
    {
    return 1;
    }
OTB Bot's avatar
STYLE    
OTB Bot committed
980
  return n * factorial(n - 1);
981
}
OTB Bot's avatar
STYLE    
OTB Bot committed
982

983
984
/** Convenience internal method */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
985
986
long ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::doubleFactorial(long n)
987
{
OTB Bot's avatar
STYLE    
OTB Bot committed
988
  if (n < 2)
989
990
991
    {
    return 1;
    }
OTB Bot's avatar
STYLE    
OTB Bot committed
992
  return n * doubleFactorial(n - 2);
993
}
OTB Bot's avatar
STYLE    
OTB Bot committed
994

995
996
/** Convenience internal method  */
template <class TLabelObject, class TLabelImage>
OTB Bot's avatar
STYLE    
OTB Bot committed
997
998
double ShapeAttributesLabelObjectFunctor<TLabelObject, TLabelImage>
::gammaN2p1(long n)
999
{
1000
  bool even = n % 2 == 0;
For faster browsing, not all history is shown. View entire blame