File: itkImageMomentsTest.cxx

package info (click to toggle)
insighttoolkit5 5.4.3-5
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 704,384 kB
  • sloc: cpp: 783,592; ansic: 628,724; xml: 44,704; fortran: 34,250; python: 22,874; sh: 4,078; pascal: 2,636; lisp: 2,158; makefile: 464; yacc: 328; asm: 205; perl: 203; lex: 146; tcl: 132; javascript: 98; csh: 81
file content (275 lines) | stat: -rw-r--r-- 9,453 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#include "itkImageMomentsCalculator.h"

#include "itkImageMaskSpatialObject.h"
#include "itkTestingMacros.h"

using PixelType = unsigned short;
using VectorType = itk::Vector<double, 3>;
using MatrixType = itk::Matrix<double, 3>;
using ImageType = itk::Image<PixelType, 3>;
using CalculatorType = itk::ImageMomentsCalculator<ImageType>;
using AffineTransformType = CalculatorType::AffineTransformType;


int
itkImageMomentsTest(int argc, char * argv[])
{
  if (argc != 2)
  {
    std::cerr << "Missing parameters." << std::endl;
    std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv);
    std::cerr << " <mask|nomask>" << std::endl;
    return EXIT_FAILURE;
  }
  const std::string maskCondition{ argv[1] };

  /* Define acceptable (absolute) error in computed results.
     All the calculations are done in double and are well-conditioned,
     so we should be able to get within a few epsilon of the right
     values.  So choose maxerr to be 10*epsilon for IEEE 754 double.
     FIXME: For some reason as yet undetermined, the Intel compiler
     produces results that are off by 12*epsilon.  This is still
     reasonably close but might deserve investigation some day when all
     the worse problems have been fixed. */
  //    double maxerr = 1.9e-15;
  constexpr double maxerr = 5.0e-15;

  /* Define the image size and physical coordinates */
  itk::Size<3> size = { { 20, 40, 80 } };
  double       origin[3] = { 0.5, 0.5, 0.5 };
  double       spacing[3] = { 0.1, 0.05, 0.025 };

  /* Define positions of the test masses in index coordinates */
  unsigned short                mass = 1; // Test mass
  itk::Index<3>::IndexValueType point[8][3] = {
    { 10 + 8, 20 + 12, 40 + 0 }, { 10 - 8, 20 - 12, 40 - 0 }, { 10 + 3, 20 - 8, 40 + 0 },
    { 10 - 3, 20 + 8, 40 - 0 },  { 10 + 0, 20 + 0, 40 + 10 }, { 10 - 0, 20 - 0, 40 - 10 },
  };

  /* Define the expected (true) results for comparison */
  double ttm = 6.0; // Total mass
  double pad[3][3] = {
    // Principal axes
    { 0.0, 0.0, 1.0 },
    { 0.6, -0.8, 0.0 },
    { 0.8, 0.6, 0.0 },
  };

  VectorType tcg;
  tcg[0] = 1.5;
  tcg[1] = 1.5;
  tcg[2] = 1.5;

  VectorType tpm;
  tpm[0] = 0.125;
  tpm[1] = 0.5;
  tpm[2] = 2.0; // Principal moments

  MatrixType tpa;
  tpa.GetVnlMatrix().set((double *)pad);

  /* Allocate a simple test image */
  auto image = ImageType::New();

  ImageType::RegionType region;
  region.SetSize(size);
  image->SetRegions(region);

  /* Set origin and spacing of physical coordinates */
  image->SetOrigin(origin);
  image->SetSpacing(spacing);
  image->Allocate();

  image->FillBuffer(PixelType{});

  /* Set a few mass points within the image */
  /* FIXME: The method used here to set the points is klutzy,
     but appears to be the only method currently supported. */
  itk::Index<3> index; /* Index over pixels */
  for (int i = 0; i < 6; ++i)
  {
    index.SetIndex(point[i]);
    image->SetPixel(index, mass);
  }

  /* Compute the moments */
  auto moments = CalculatorType::New();
  moments->SetImage(image);
  if (maskCondition == std::string("mask"))
  {
    // Test the mask spatial object for masked ImageMomentsTest
    // Make a mask that covers the entire image space
    using MaskImageType = itk::Image<unsigned char, 3>;
    auto maskimg = MaskImageType::New();
    maskimg->CopyInformation(image);
    maskimg->SetRegions(image->GetLargestPossibleRegion());
    maskimg->Allocate();
    // Masking the entire image should not change the computation results.
    maskimg->FillBuffer(itk::NumericTraits<PixelType>::OneValue());

    // convert mask image to mask
    using LFFImageMaskSpatialObjectType = typename itk::ImageMaskSpatialObject<MaskImageType::ImageDimension>;
    auto mask = LFFImageMaskSpatialObjectType::New();
    mask->SetImage(maskimg.GetPointer());
    mask->Update();
    // Purposefully use the base class type
    typename itk::SpatialObject<MaskImageType::ImageDimension>::Pointer test =
      dynamic_cast<itk::SpatialObject<MaskImageType::ImageDimension> *>(mask.GetPointer());
    if (test.IsNull())
    {
      itkGenericExceptionMacro("Failed conversion to SpatialObject base class.");
    }
    moments->SetSpatialObjectMask(test.GetPointer());
  }
  moments->Compute();

  /* Printout info */
  moments->Print(std::cout);

  double     ctm = moments->GetTotalMass();
  VectorType ccg = moments->GetCenterOfGravity();
  VectorType cpm = moments->GetPrincipalMoments();
  MatrixType cpa = moments->GetPrincipalAxes();

  /* Flip the principal axes if necessary.

  The eigenvector solution is still valid if multiplied by a constant.
  Since the eigenvectors (principal axes) are normalized, this constant
  can be -1.  */
  if (cpa(1, 0) < 0.0) // Should be 0.6
  {
    for (unsigned int row = 0; row < 3; ++row)
    {
      cpa(row, 0) *= -1;
    }
  }
  if (cpa(1, 1) > 0.0) // Should be -0.8
  {
    for (unsigned int row = 0; row < 3; ++row)
    {
      cpa(row, 1) *= -1;
    }
  }

  /* Report the various non-central moments */
  // FIXME:  Indentation is not handled correctly in matrix output
  std::cout << "\nTotal mass = " << ctm << std::endl;
  std::cout << "True total mass = " << ttm << std::endl;
  std::cout << "\nFirst moments about index origin =\n";
  std::cout << "   " << moments->GetFirstMoments() << std::endl;
  std::cout << "\nSecond moments about index origin =\n";
  std::cout << "   " << moments->GetSecondMoments() << std::endl;

  /* Report the center of gravity and central moments */
  std::cout << "\nCenter of gravity =\n";
  std::cout << "   " << ccg << '\n';
  std::cout << "True center of gravity =\n";
  std::cout << "   " << tcg << '\n';
  std::cout << "\nSecond central moments =\n";
  std::cout << "   " << moments->GetCentralMoments() << '\n';

  /* Report principal moments and axes */
  std::cout << "\nPrincipal moments = \n";
  std::cout << "   " << cpm << '\n';
  std::cout << "True principal moments = \n";
  std::cout << "   " << tpm << '\n';
  std::cout << "\nPrincipal axes = \n";
  std::cout << "   " << cpa << '\n';
  std::cout << "True principal axes = \n";
  std::cout << "   " << tpa << '\n';

  /* Compute transforms between principal and physical axes */
  /* FIXME: Automatically check correctness of these results? */
  AffineTransformType::Pointer pa2p = moments->GetPrincipalAxesToPhysicalAxesTransform();
  std::cout << "\nPrincipal axes to physical axes transform:\n";
  std::cout << pa2p->GetMatrix() << std::endl;
  AffineTransformType::Pointer p2pa = moments->GetPhysicalAxesToPrincipalAxesTransform();
  std::cout << "\nPhysical axes to principal axes transform:\n";
  std::cout << p2pa->GetMatrix() << std::endl;

  /* Do some error checking on the transforms */
  double dist = pa2p->Metric(pa2p);
  std::cout << "Distance from self to self = " << dist << std::endl;
  auto p2pa2p = AffineTransformType::New();
  p2pa2p->Compose(p2pa);
  p2pa2p->Compose(pa2p);
  double trerr = p2pa2p->Metric();
  std::cout << "Distance from composition to identity = ";
  std::cout << trerr << std::endl;


  /* Compute and report max abs error in computed */
  double tmerr = itk::Math::abs(ttm - ctm); // Error in total mass
  double cgerr = 0.0;                       // Error in center of gravity
  double pmerr = 0.0;                       // Error in moments
  double paerr = 0.0;                       // Error in axes

  for (int i = 0; i < 3; ++i)
  {
    if (itk::Math::abs(ccg[i] - tcg[i]) > cgerr)
    {
      cgerr = itk::Math::abs(ccg[i] - tcg[i]);
    }
    if (itk::Math::abs(cpm[i] - tpm[i]) > pmerr)
    {
      pmerr = itk::Math::abs(cpm[i] - tpm[i]);
    }
    for (int j = 0; j < 3; ++j)
    {
      if (itk::Math::abs(cpa[i][j] - tpa[i][j]) > paerr)
      {
        paerr = itk::Math::abs(cpa[i][j] - tpa[i][j]);
      }
    }
  }

  std::cout << "\nErrors found in:\n";
  std::cout << "   Total mass        = " << tmerr << std::endl;
  std::cout << "   Center of gravity = " << cgerr << std::endl;
  std::cout << "   Principal moments = " << pmerr << std::endl;
  std::cout << "   Principal axes    = " << paerr << std::endl;
  std::cout << "   Transformations   = " << trerr << std::endl;

  /* Return error if differences are too large */
  int stat = tmerr > maxerr || cgerr > maxerr || pmerr > maxerr || paerr > maxerr || trerr > maxerr;

  std::cout << std::endl;
  bool pass;
  if (stat)
  {
    std::cout << "Errors are larger than defined maximum value." << std::endl;
    std::cout << "Test FAILED !" << std::endl;
    pass = false;
  }
  else
  {
    std::cout << "Errors are acceptable" << std::endl;
    std::cout << "Test PASSED !" << std::endl;
    pass = true;
  }

  if (!pass)
  {
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}