| .. default-domain:: cpp |
| |
| .. cpp:namespace:: ceres |
| |
| .. _chapter-interfacing_with_automatic_differentiation: |
| |
| Interfacing with Automatic Differentiation |
| ========================================== |
| |
| Automatic differentiation is straightforward to use in cases where an |
| explicit expression for the cost function is available. But this is |
| not always possible. Often one has to interface with external routines |
| or data. In this chapter we will consider a number of different ways |
| of doing so. |
| |
| To do this, we will consider the problem of finding parameters |
| :math:`\theta` and :math:`t` that solve an optimization problem of the |
| form: |
| |
| .. math:: |
| \min & \quad \sum_i \left \|y_i - f\left (\|q_{i}\|^2\right) q_i |
| \right \|^2\\ |
| \text{such that} & \quad q_i = R(\theta) x_i + t |
| |
| Here, :math:`R` is a two dimensional rotation matrix parameterized |
| using the angle :math:`\theta` and :math:`t` is a two dimensional |
| vector. :math:`f` is an external distortion function. |
| |
| We begin by considering the case, where we have a templated function |
| :code:`TemplatedComputeDistortion` that can compute the function |
| :math:`f`. Then the implementation of the corresponding residual |
| functor is straightforward and will look as follows: |
| |
| .. code-block:: c++ |
| :emphasize-lines: 21 |
| |
| template <typename T> T TemplatedComputeDistortion(const T r2) { |
| const double k1 = 0.0082; |
| const double k2 = 0.000023; |
| return 1.0 + k1 * y2 + k2 * r2 * r2; |
| } |
| |
| struct Affine2DWithDistortion { |
| Affine2DWithDistortion(const double x_in[2], const double y_in[2]) { |
| x[0] = x_in[0]; |
| x[1] = x_in[1]; |
| y[0] = y_in[0]; |
| y[1] = y_in[1]; |
| } |
| |
| template <typename T> |
| bool operator()(const T* theta, |
| const T* t, |
| T* residuals) const { |
| const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0]; |
| const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1]; |
| const T f = TemplatedComputeDistortion(q_0 * q_0 + q_1 * q_1); |
| residuals[0] = y[0] - f * q_0; |
| residuals[1] = y[1] - f * q_1; |
| return true; |
| } |
| |
| double x[2]; |
| double y[2]; |
| }; |
| |
| So far so good, but let us now consider three ways of defining |
| :math:`f` which are not directly amenable to being used with automatic |
| differentiation: |
| |
| #. A non-templated function that evaluates its value. |
| #. A function that evaluates its value and derivative. |
| #. A function that is defined as a table of values to be interpolated. |
| |
| We will consider them in turn below. |
| |
| A function that returns its value |
| ---------------------------------- |
| |
| Suppose we were given a function :code:`ComputeDistortionValue` with |
| the following signature |
| |
| .. code-block:: c++ |
| |
| double ComputeDistortionValue(double r2); |
| |
| that computes the value of :math:`f`. The actual implementation of the |
| function does not matter. Interfacing this function with |
| :code:`Affine2DWithDistortion` is a three step process: |
| |
| 1. Wrap :code:`ComputeDistortionValue` into a functor |
| :code:`ComputeDistortionValueFunctor`. |
| 2. Numerically differentiate :code:`ComputeDistortionValueFunctor` |
| using :class:`NumericDiffCostFunction` to create a |
| :class:`CostFunction`. |
| 3. Wrap the resulting :class:`CostFunction` object using |
| :class:`CostFunctionToFunctor`. The resulting object is a functor |
| with a templated :code:`operator()` method, which pipes the |
| Jacobian computed by :class:`NumericDiffCostFunction` into the |
| appropriate :code:`Jet` objects. |
| |
| An implementation of the above three steps looks as follows: |
| |
| .. code-block:: c++ |
| :emphasize-lines: 15,16,17,18,19,20, 29 |
| |
| struct ComputeDistortionValueFunctor { |
| bool operator()(const double* r2, double* value) const { |
| *value = ComputeDistortionValue(r2[0]); |
| return true; |
| } |
| }; |
| |
| struct Affine2DWithDistortion { |
| Affine2DWithDistortion(const double x_in[2], const double y_in[2]) { |
| x[0] = x_in[0]; |
| x[1] = x_in[1]; |
| y[0] = y_in[0]; |
| y[1] = y_in[1]; |
| |
| compute_distortion.reset(new ceres::CostFunctionToFunctor<1, 1>( |
| new ceres::NumericDiffCostFunction<ComputeDistortionValueFunctor, |
| ceres::CENTRAL, |
| 1, |
| 1>( |
| new ComputeDistortionValueFunctor))); |
| } |
| |
| template <typename T> |
| bool operator()(const T* theta, const T* t, T* residuals) const { |
| const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0]; |
| const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1]; |
| const T r2 = q_0 * q_0 + q_1 * q_1; |
| T f; |
| (*compute_distortion)(&r2, &f); |
| residuals[0] = y[0] - f * q_0; |
| residuals[1] = y[1] - f * q_1; |
| return true; |
| } |
| |
| double x[2]; |
| double y[2]; |
| std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion; |
| }; |
| |
| |
| A function that returns its value and derivative |
| ------------------------------------------------ |
| |
| Now suppose we are given a function :code:`ComputeDistortionValue` |
| thatis able to compute its value and optionally its Jacobian on demand |
| and has the following signature: |
| |
| .. code-block:: c++ |
| |
| void ComputeDistortionValueAndJacobian(double r2, |
| double* value, |
| double* jacobian); |
| |
| Again, the actual implementation of the function does not |
| matter. Interfacing this function with :code:`Affine2DWithDistortion` |
| is a two step process: |
| |
| 1. Wrap :code:`ComputeDistortionValueAndJacobian` into a |
| :class:`CostFunction` object which we call |
| :code:`ComputeDistortionFunction`. |
| 2. Wrap the resulting :class:`ComputeDistortionFunction` object using |
| :class:`CostFunctionToFunctor`. The resulting object is a functor |
| with a templated :code:`operator()` method, which pipes the |
| Jacobian computed by :class:`NumericDiffCostFunction` into the |
| appropriate :code:`Jet` objects. |
| |
| The resulting code will look as follows: |
| |
| .. code-block:: c++ |
| :emphasize-lines: 21,22, 33 |
| |
| class ComputeDistortionFunction : public ceres::SizedCostFunction<1, 1> { |
| public: |
| virtual bool Evaluate(double const* const* parameters, |
| double* residuals, |
| double** jacobians) const { |
| if (!jacobians) { |
| ComputeDistortionValueAndJacobian(parameters[0][0], residuals, NULL); |
| } else { |
| ComputeDistortionValueAndJacobian(parameters[0][0], residuals, jacobians[0]); |
| } |
| return true; |
| } |
| }; |
| |
| struct Affine2DWithDistortion { |
| Affine2DWithDistortion(const double x_in[2], const double y_in[2]) { |
| x[0] = x_in[0]; |
| x[1] = x_in[1]; |
| y[0] = y_in[0]; |
| y[1] = y_in[1]; |
| compute_distortion.reset( |
| new ceres::CostFunctionToFunctor<1, 1>(new ComputeDistortionFunction)); |
| } |
| |
| template <typename T> |
| bool operator()(const T* theta, |
| const T* t, |
| T* residuals) const { |
| const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0]; |
| const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1]; |
| const T r2 = q_0 * q_0 + q_1 * q_1; |
| T f; |
| (*compute_distortion)(&r2, &f); |
| residuals[0] = y[0] - f * q_0; |
| residuals[1] = y[1] - f * q_1; |
| return true; |
| } |
| |
| double x[2]; |
| double y[2]; |
| std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion; |
| }; |
| |
| |
| A function that is defined as a table of values |
| ----------------------------------------------- |
| |
| The third and final case we will consider is where the function |
| :math:`f` is defined as a table of values on the interval :math:`[0, |
| 100)`, with a value for each integer. |
| |
| .. code-block:: c++ |
| |
| vector<double> distortion_values; |
| |
| There are many ways of interpolating a table of values. Perhaps the |
| simplest and most common method is linear interpolation. But it is not |
| a great idea to use linear interpolation because the interpolating |
| function is not differentiable at the sample points. |
| |
| A simple (well behaved) differentiable interpolation is the `Cubic |
| Hermite Spline |
| <http://en.wikipedia.org/wiki/Cubic_Hermite_spline>`_. Ceres Solver |
| ships with routines to perform Cubic & Bi-Cubic interpolation that is |
| automatic differentiation friendly. |
| |
| Using Cubic interpolation requires first constructing a |
| :class:`Grid1D` object to wrap the table of values and then |
| constructing a :class:`CubicInterpolator` object using it. |
| |
| The resulting code will look as follows: |
| |
| .. code-block:: c++ |
| :emphasize-lines: 10,11,12,13, 24, 32,33 |
| |
| struct Affine2DWithDistortion { |
| Affine2DWithDistortion(const double x_in[2], |
| const double y_in[2], |
| const std::vector<double>& distortion_values) { |
| x[0] = x_in[0]; |
| x[1] = x_in[1]; |
| y[0] = y_in[0]; |
| y[1] = y_in[1]; |
| |
| grid.reset(new ceres::Grid1D<double, 1>( |
| &distortion_values[0], 0, distortion_values.size())); |
| compute_distortion.reset( |
| new ceres::CubicInterpolator<ceres::Grid1D<double, 1> >(*grid)); |
| } |
| |
| template <typename T> |
| bool operator()(const T* theta, |
| const T* t, |
| T* residuals) const { |
| const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0]; |
| const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1]; |
| const T r2 = q_0 * q_0 + q_1 * q_1; |
| T f; |
| compute_distortion->Evaluate(r2, &f); |
| residuals[0] = y[0] - f * q_0; |
| residuals[1] = y[1] - f * q_1; |
| return true; |
| } |
| |
| double x[2]; |
| double y[2]; |
| std::unique_ptr<ceres::Grid1D<double, 1> > grid; |
| std::unique_ptr<ceres::CubicInterpolator<ceres::Grid1D<double, 1> > > compute_distortion; |
| }; |
| |
| In the above example we used :class:`Grid1D` and |
| :class:`CubicInterpolator` to interpolate a one dimensional table of |
| values. :class:`Grid2D` combined with :class:`CubicInterpolator` lets |
| the user to interpolate two dimensional tables of values. Note that |
| neither :class:`Grid1D` or :class:`Grid2D` are limited to scalar |
| valued functions, they also work with vector valued functions. |