← Back to branch summary

~ubuntu-branches/ubuntu/utopic/blender/utopic-proposed

« back to all changes in this revision

Viewing changes to extern/libmv/third_party/ceres/internal/ceres/corrector.cc

  • Committer: Package Import Robot
  • Author(s): Matthias Klose
  • Date: 2014-02-19 11:24:23 UTC
  • mfrom: (14.2.23 sid)
  • Revision ID: package-import@ubuntu.com-20140219112423-rkmaz2m7ha06d4tk
Tags: 2.69-3ubuntu1
* Merge with Debian; remaining changes:
  - Configure without OpenImageIO on armhf, as it is not available on
    Ubuntu.

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/libmv/third_party/ceres/internal/ceres/corrector.cc
32
32
 
33
33
#include <cstddef>
34
34
#include <cmath>
35
 
#include "ceres/internal/eigen.h"
36
35
#include "glog/logging.h"
37
36
 
38
37
namespace ceres {
90
89
  // 0.5 *  alpha^2 - alpha - rho'' / rho' *  z'z = 0.
91
90
  //
92
91
  // Start by calculating the discriminant D.
93
 
  const double D = 1.0 + 2.0 * sq_norm*rho[2] / rho[1];
 
92
  const double D = 1.0 + 2.0 * sq_norm * rho[2] / rho[1];
94
93
 
95
94
  // Since both rho[1] and rho[2] are guaranteed to be positive at
96
95
  // this point, we know that D > 1.0.
102
101
  alpha_sq_norm_ = alpha / sq_norm;
103
102
}
104
103
 
105
 
void Corrector::CorrectResiduals(int nrow, double* residuals) {
 
104
void Corrector::CorrectResiduals(int num_rows, double* residuals) {
106
105
  DCHECK(residuals != NULL);
107
 
  VectorRef r_ref(residuals, nrow);
108
106
  // Equation 11 in BANS.
109
 
  r_ref *= residual_scaling_;
 
107
  for (int r = 0; r < num_rows; ++r) {
 
108
    residuals[r] *= residual_scaling_;
 
109
  }
110
110
}
111
111
 
112
 
void Corrector::CorrectJacobian(int nrow, int ncol,
113
 
                                double* residuals, double* jacobian) {
 
112
void Corrector::CorrectJacobian(int num_rows,
 
113
                                int num_cols,
 
114
                                double* residuals,
 
115
                                double* jacobian) {
114
116
  DCHECK(residuals != NULL);
115
117
  DCHECK(jacobian != NULL);
116
 
 
117
 
  if (nrow == 1) {
118
 
    // Specialization for the case where the residual is a scalar.
119
 
    VectorRef j_ref(jacobian, ncol);
120
 
    j_ref *= sqrt_rho1_ * (1.0 - alpha_sq_norm_ * pow(*residuals, 2));
121
 
  } else {
122
 
    ConstVectorRef r_ref(residuals, nrow);
123
 
    MatrixRef j_ref(jacobian, nrow, ncol);
124
 
 
125
 
    // Equation 11 in BANS.
126
 
    j_ref = sqrt_rho1_ * (j_ref - alpha_sq_norm_ *
127
 
                          r_ref * (r_ref.transpose() * j_ref));
 
118
  // Equation 11 in BANS.
 
119
  //
 
120
  //  J = sqrt(rho) * (J - alpha^2 r * r' J)
 
121
  //
 
122
  // In days gone by this loop used to be a single Eigen expression of
 
123
  // the form
 
124
  //
 
125
  //  J = sqrt_rho1_ * (J - alpha_sq_norm_ * r* (r.transpose() * J));
 
126
  //
 
127
  // Which turns out to about 17x slower on bal problems. The reason
 
128
  // is that Eigen is unable to figure out that this expression can be
 
129
  // evaluated columnwise and ends up creating a temporary.
 
130
  for (int c = 0; c < num_cols; ++c) {
 
131
    double r_transpose_j = 0.0;
 
132
    for (int r = 0; r < num_rows; ++r) {
 
133
      r_transpose_j += jacobian[r * num_cols + c] * residuals[r];
 
134
    }
 
135
 
 
136
    for (int r = 0; r < num_rows; ++r) {
 
137
      jacobian[r * num_cols + c] = sqrt_rho1_ *
 
138
          (jacobian[r * num_cols + c] -
 
139
           alpha_sq_norm_ * residuals[r] * r_transpose_j);
 
140
    }
128
141
  }
129
142
}
130
143