Seg3D  2.4
Seg3D is a free volume segmentation and processing tool developed by the NIH Center for Integrative Biomedical Computing at the University of Utah Scientific Computing and Imaging (SCI) Institute.
itkNumericInverse.h
1 /*
2  For more information, please see: http://software.sci.utah.edu
3 
4  The MIT License
5 
6  Copyright (c) 2016 Scientific Computing and Imaging Institute,
7  University of Utah.
8 
9 
10  Permission is hereby granted, free of charge, to any person obtaining a
11  copy of this software and associated documentation files (the "Software"),
12  to deal in the Software without restriction, including without limitation
13  the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  and/or sell copies of the Software, and to permit persons to whom the
15  Software is furnished to do so, subject to the following conditions:
16 
17  The above copyright notice and this permission notice shall be included
18  in all copies or substantial portions of the Software.
19 
20  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  DEALINGS IN THE SOFTWARE.
27 */
28 
29 // File : itkNumericTransform.h
30 // Author : Pavel A. Koshevoy
31 // Created : 2005/06/03 10:16
32 // Copyright : (C) 2004-2008 University of Utah
33 // Description : A numerical inverse transform class, based on the
34 // Newton-Raphson method for nonlinear systems of equations.
35 
36 #ifndef __itkNumericInverse_h
37 #define __itkNumericInverse_h
38 
39 // system includes:
40 #include <iostream>
41 #include <vector>
42 
43 // ITK includes:
44 #include <itkTransform.h>
45 #include <itkMacro.h>
46 
47 // VXL includes:
48 #include <vnl/algo/vnl_svd.h>
49 
50 namespace help
51 {
52 
53 //----------------------------------------------------------------
54 // nonlinear_system_evaluator_t
55 //
56 template <class ScalarType>
58 {
59 public:
60  virtual ~nonlinear_system_evaluator_t() {}
61 
62  // evaluate the nonlinear system of equations
63  // and its Jacobian (dF/dx) at a given point:
64  virtual void
65  operator() (const std::vector<ScalarType> & x,
66  std::vector<ScalarType> & F,
67  std::vector<std::vector<ScalarType> > & J) const = 0;
68 };
69 
70 //----------------------------------------------------------------
71 // NewtonRaphson
72 //
73 template <class ScalarType>
74 bool
75 NewtonRaphson(const nonlinear_system_evaluator_t<ScalarType> & usrfun,
76  std::vector<ScalarType> & x, // estimated root point
77  const unsigned int & ntrial, // number of iterations
78  const ScalarType & tolx, // convergence tolerance in x
79  const ScalarType & tolf) // convergence tolerance in F
80 {
81  // shortcut:
82  const unsigned int n = x.size();
83 
84  std::vector<ScalarType> F(n);
85  std::vector<std::vector<ScalarType> > J(n);
86  for (unsigned int i = 0; i < n; i++) J[i].resize(n);
87 
88  vnl_matrix<ScalarType> A(n, n);
89  vnl_vector<ScalarType> b(n);
90 
91  for (unsigned int k = 0; k < ntrial; k++)
92  {
93  // evaluate the function at the current position:
94  usrfun(x, F, J);
95 
96  // check for function convergence:
97  ScalarType errf = ScalarType(0);
98  for (unsigned int i = 0; i < n; i++)
99  {
100  errf += fabs(F[i]);
101  }
102  if (errf <= tolf) break;
103 
104  // setup the left hand side of the linear system:
105  for (unsigned int i = 0; i < n; i++)
106  {
107  for (unsigned int j = 0; j < n; j++)
108  {
109  A[i][j] = J[i][j];
110  }
111  }
112 
113  // setup the right hand side of the linear system:
114  for (unsigned int i = 0; i < n; i++)
115  {
116  b[i] = -F[i];
117  }
118 
119  vnl_svd<double> svd(A);
120  vnl_vector<double> dx = svd.solve(b);
121 
122  // check for root convergence:
123  ScalarType errx = ScalarType(0);
124  for (unsigned int i = 0; i < n; i++)
125  {
126  errx += fabs(dx[i]);
127  x[i] += dx[i];
128  }
129  if (errx <= tolx) break;
130  }
131 
132  return true;
133 }
134 
135 } // namespace help
136 
137 namespace itk
138 {
139 
140 
141 //----------------------------------------------------------------
142 // NumericInverse
143 //
144 template <class TTransform>
146  public help::nonlinear_system_evaluator_t<typename TTransform::ScalarType>
147 {
148 public:
149  typedef TTransform TransformType;
150  typedef typename TTransform::ScalarType ScalarType;
151 
152  NumericInverse(const TransformType & transform):
153  transform_(transform)
154  {}
155 
156  // virtual:
157  virtual
158  void operator() (const std::vector<ScalarType> & x,
159  std::vector<ScalarType> & F,
160  std::vector<std::vector<ScalarType> > & J) const
161  {
162  transform_.eval(x, F, J);
163 
164  const unsigned int n = x.size();
165  for (unsigned int i = 0; i < n; i++)
166  {
167  F[i] -= y_[i];
168  }
169  }
170 
171  // If y = Transform(x), then x = BackTransform(y).
172  // given y, find x:
173  bool transform(const std::vector<ScalarType> & y,
174  std::vector<ScalarType> & x,
175  bool x_is_initialized = false) const
176  {
177  y_ = y;
178  if (!x_is_initialized) x = y;
179  return NewtonRaphson(*this, x, 50, 1e-12, 1e-12);
180  }
181 
182 private:
183  // the transform whose inverse we are trying to evaluate:
184  const TransformType & transform_;
185 
186  // the point for which we are tryying to find the inverse mapping:
187  mutable std::vector<ScalarType> y_;
188 };
189 
190 } // namespace itk
191 
192 #endif // __itkNumericInverse_h
Definition: itkNormalizeImageFilterWithMask.h:48
Definition: itkNumericInverse.h:57
Definition: itkNumericInverse.h:145
Definition: itkNumericInverse.h:50