1 #ifndef rgrl_est_proj_func_txx_
2 #define rgrl_est_proj_func_txx_
3 //:
4 // \file
5 // \author Gehua Yang
6 // \date March 2007
7 // \brief a generic class to estimate a homogeneous projection matrix using L-M
8
9 #include "rgrl_est_proj_func.h"
10
11 #include <rgrl/rgrl_estimator.h>
12 #include <rgrl/rgrl_match_set.h>
13 #include <rgrl/rgrl_set_of.h>
14
15 #include <vnl/vnl_matrix_fixed.h>
16 #include <vnl/vnl_vector_fixed.h>
17 #include <vnl/vnl_math.h>
18 #include <vnl/vnl_numeric_traits.h>
19 #include <vnl/vnl_transpose.h>
20 #include <vnl/algo/vnl_levenberg_marquardt.h>
21 #include <vnl/algo/vnl_svd.h>
22
23 #include <cassert>
24 #ifdef _MSC_VER
25 # include <vcl_msvc_warnings.h>
26 #endif
27
28 const static unsigned int maxval_unsigned = vnl_numeric_traits<unsigned int>::maxval;
29
30 template <unsigned int Tdim, unsigned int Fdim>
31 rgrl_est_proj_func<Tdim, Fdim>::
rgrl_est_proj_func(rgrl_set_of<rgrl_match_set_sptr> const & matches,bool with_grad)32 rgrl_est_proj_func( rgrl_set_of<rgrl_match_set_sptr> const& matches,
33 bool with_grad )
34 : vnl_least_squares_function( (Fdim+1)*(Tdim+1)-1,
35 rgrl_est_matches_residual_number(matches),
36 with_grad ? use_gradient : no_gradient ),
37 matches_ptr_( &matches ),
38 from_centre_(double(0)), to_centre_(double(0)),
39 index_row_(maxval_unsigned), index_col_(maxval_unsigned),
40 max_num_iterations_(50),
41 relative_threshold_(1e-7),
42 zero_svd_thres_(1e-5)
43 {
44 }
45
46 template <unsigned int Tdim, unsigned int Fdim>
rgrl_est_proj_func(bool with_grad)47 rgrl_est_proj_func<Tdim, Fdim>::rgrl_est_proj_func(bool with_grad)
48 : vnl_least_squares_function((Fdim + 1) * (Tdim + 1) - 1,
49 1000 /*artificial number to avoid warning*/,
50 with_grad ? use_gradient : no_gradient),
51
52 from_centre_(double(0)), to_centre_(double(0)),
53 index_row_(maxval_unsigned), index_col_(maxval_unsigned)
54
55 {}
56
57 //: convert parameters
58 template <unsigned int Tdim, unsigned int Fdim>
59 void
60 rgrl_est_proj_func<Tdim, Fdim>::
convert_parameters(vnl_vector<double> & params,vnl_matrix_fixed<double,Tdim+1,Fdim+1> proj_matrix,vnl_vector_fixed<double,Fdim> const & fc,vnl_vector_fixed<double,Tdim> const & tc)61 convert_parameters( vnl_vector<double>& params,
62 vnl_matrix_fixed<double, Tdim+1, Fdim+1> proj_matrix,
63 vnl_vector_fixed<double, Fdim> const& fc,
64 vnl_vector_fixed<double, Tdim> const& tc )
65 {
66 // set centres
67 from_centre_ = fc;
68 to_centre_ = tc;
69
70 // make the prejection matrix to centre around from_centre_
71 vnl_matrix_fixed<double, Fdim+1, Fdim+1> from_centre_matrix;
72 from_centre_matrix.set_identity();
73 for ( unsigned i=0; i<Fdim; ++i )
74 from_centre_matrix( i, Fdim ) = from_centre_[i];
75
76 vnl_matrix_fixed<double, Tdim+1, Tdim+1> to_centre_matrix;
77 to_centre_matrix.set_identity();
78 for ( unsigned i=0; i<Tdim; ++i )
79 to_centre_matrix( i, Tdim ) = -to_centre_[i];
80
81 proj_matrix = to_centre_matrix * proj_matrix * from_centre_matrix;
82
83 // find out which element in the projection matrix has largest value
84 double max_val = 0;
85 index_row_ = index_col_ = -1; // init to bad value
86 for ( unsigned i=0; i<Tdim+1; ++i )
87 for ( unsigned j=0; j<Fdim+1; ++j )
88 if ( std::abs( proj_matrix(i,j) ) > max_val ) {
89 index_row_ = i;
90 index_col_ = j;
91 max_val = std::abs( proj_matrix(i,j) );
92 }
93
94 // normalize the proj_matrix to have the largest value as 1
95 proj_matrix /= proj_matrix( index_row_, index_col_ );
96
97 // fill in params
98 params.set_size( proj_size_-1 );
99 for ( unsigned k=0,i=0; i<Tdim+1; ++i )
100 for ( unsigned j=0; j<Fdim+1; ++j ) {
101 if ( i==index_row_ && j==index_col_ ) {
102 continue;
103 }
104 // fill in elements in order
105 params[k++] = proj_matrix(i,j);
106 }
107 }
108
109 //: uncenter projection matrix
110 template <unsigned int Tdim, unsigned int Fdim>
111 vnl_matrix_fixed<double, Tdim+1, Fdim+1>
112 rgrl_est_proj_func<Tdim, Fdim>::
uncentre_proj(vnl_matrix_fixed<double,Tdim+1,Fdim+1> const & proj) const113 uncentre_proj( vnl_matrix_fixed<double, Tdim+1, Fdim+1> const& proj ) const
114 {
115 // make the prejection matrix to centre around from_centre_
116 vnl_matrix_fixed<double, Fdim+1, Fdim+1> from_centre_matrix;
117 from_centre_matrix.set_identity();
118 for ( unsigned i=0; i<Fdim; ++i )
119 from_centre_matrix( i, Fdim ) = - from_centre_[i];
120
121 vnl_matrix_fixed<double, Tdim+1, Tdim+1> to_centre_matrix;
122 to_centre_matrix.set_identity();
123 for ( unsigned i=0; i<Tdim; ++i )
124 to_centre_matrix( i, Tdim ) = to_centre_[i];
125
126 return to_centre_matrix * proj * from_centre_matrix;
127 }
128
129 template <unsigned int Tdim, unsigned int Fdim>
130 void
131 rgrl_est_proj_func<Tdim, Fdim>::
full_proj_jacobian(vnl_matrix_fixed<double,Tdim,(Fdim+1)* (Tdim+1)> & complete_jac,vnl_matrix_fixed<double,Tdim+1,Fdim+1> const & proj,vnl_vector_fixed<double,Fdim> const & from) const132 full_proj_jacobian( vnl_matrix_fixed<double, Tdim, (Fdim+1)*(Tdim+1)>& complete_jac,
133 vnl_matrix_fixed<double, Tdim+1, Fdim+1> const& proj,
134 vnl_vector_fixed<double, Fdim> const& from ) const
135 {
136 // subtract centre
137 const vnl_vector_fixed<double, Fdim> from_centred = from-from_centre_;
138
139 vnl_vector_fixed<double, Tdim+1> homo;
140 vnl_matrix_fixed<double, Tdim+1, proj_size_> jf(0.0); // grad in homogeneous coordinate
141 vnl_matrix_fixed<double, Tdim, Tdim+1> jg(0.0); // grad of division, [u/w, v/w]^T
142
143 rgrl_est_proj_map_homo_point<Tdim, Fdim>( homo, proj, from_centred );
144
145 // 1. linear gradient in homogeneous coordinate
146 // fill jf (linear gradient) with 1.0 on elements corresponding to shift
147 for ( unsigned i=0; i<Tdim+1; ++i )
148 jf( i, i*(Fdim+1)+Fdim ) = 1.0;
149
150 // skip the ones corresponding to shift
151 for ( unsigned index=0,i=0; i<Tdim+1; ++i,index+=(Fdim+1) )
152 for ( unsigned j=0; j<Fdim; ++j ) {
153 jf( i, index+j ) = from_centred[j];
154 }
155
156 // 2. gradient of making division
157 const double homo_last_neg_sqr = -vnl_math::sqr(homo[Tdim]);
158 const double homo_last_div = 1.0/homo[Tdim];
159 for ( unsigned i=0; i<Tdim; ++i )
160 jg(i,i) = homo_last_div;
161 for ( unsigned i=0; i<Tdim; ++i )
162 jg(i,Tdim) = homo[i] / homo_last_neg_sqr;
163
164 // 3. complete jacobian
165 // since Jab_g(f(p)) = Jac_g * Jac_f
166 complete_jac = jg * jf;
167 }
168
169 template <unsigned int Tdim, unsigned int Fdim>
170 void
171 rgrl_est_proj_func<Tdim, Fdim>::
reduced_proj_jacobian(vnl_matrix_fixed<double,Tdim,(Fdim+1)* (Tdim+1)-1> & base_jac,vnl_matrix_fixed<double,Tdim+1,Fdim+1> const & proj,vnl_vector_fixed<double,Fdim> const & from) const172 reduced_proj_jacobian( vnl_matrix_fixed<double, Tdim, (Fdim+1)*(Tdim+1)-1>& base_jac,
173 vnl_matrix_fixed<double, Tdim+1, Fdim+1> const& proj,
174 vnl_vector_fixed<double, Fdim> const& from ) const
175 {
176 vnl_matrix_fixed<double, Tdim, proj_size_> complete_jac;
177 full_proj_jacobian( complete_jac, proj, from );
178
179 // 4. remove the element being held as constant 1
180 const unsigned index = (index_row_*(Fdim+1)+index_col_);
181 for ( unsigned i=0,j=0; i<proj_size_; ++i )
182 if ( i != index ) {
183 for ( unsigned k=0; k<Tdim; ++k )
184 base_jac(k, j) = complete_jac(k, i);
185 ++j;
186 }
187 }
188
189 template <unsigned int Tdim, unsigned int Fdim>
190 void
191 rgrl_est_proj_func<Tdim, Fdim>::
proj_jac_wrt_loc(vnl_matrix_fixed<double,Tdim,Fdim> & jac_loc,vnl_matrix_fixed<double,Tdim+1,Fdim+1> const & proj,vnl_vector_fixed<double,Fdim> const & from) const192 proj_jac_wrt_loc( vnl_matrix_fixed<double, Tdim, Fdim>& jac_loc,
193 vnl_matrix_fixed<double, Tdim+1, Fdim+1> const& proj,
194 vnl_vector_fixed<double, Fdim> const& from ) const
195 {
196 // 1. linear gradient in homogeneous coordinate
197 // fill jf (linear gradient) with 1.0 on elements corresponding to shift
198 vnl_matrix_fixed<double, Tdim+1, Fdim> jf; // grad in homogeneous coordinate
199 for ( unsigned i=0; i<Tdim+1; ++i )
200 for ( unsigned j=0; j<Fdim; ++j )
201 jf(i, j) = proj(i, j);
202
203 // 2. gradient of making division
204 vnl_vector_fixed<double, Tdim+1> homo;
205 rgrl_est_proj_map_homo_point<Tdim, Fdim>( homo, proj, from-from_centre_ );
206
207 vnl_matrix_fixed<double, Tdim, Tdim+1> jg(0.0); // grad of division, [u/w, v/w]^T
208 const double homo_last_neg_sqr = -vnl_math::sqr(homo[Tdim]);
209 const double homo_last_div = 1.0/homo[Tdim];
210 for ( unsigned i=0; i<Tdim; ++i )
211 jg(i,i) = homo_last_div;
212 for ( unsigned i=0; i<Tdim; ++i )
213 jg(i,Tdim) = homo[i] / homo_last_neg_sqr;
214
215 // 3. complete jacobian
216 // since Jab_g(f(p)) = Jac_g * Jac_f
217 jac_loc = jg * jf;
218 }
219
220 template <unsigned int Tdim, unsigned int Fdim>
221 void
222 rgrl_est_proj_func<Tdim, Fdim>::
restored_centered_proj(vnl_matrix_fixed<double,Tdim+1,Fdim+1> & proj_matrix,vnl_vector<double> const & params) const223 restored_centered_proj( vnl_matrix_fixed<double, Tdim+1, Fdim+1>& proj_matrix,
224 vnl_vector<double> const& params ) const
225 {
226 assert( params.size()+1 >= proj_size_ );
227 for ( unsigned k=0,i=0; i<Tdim+1; ++i )
228 for ( unsigned j=0; j<Fdim+1; ++j ) {
229 if ( i==index_row_ && j==index_col_ ) {
230 proj_matrix(i,j) = 1.0;
231 continue;
232 }
233 // fill in elements in order
234 proj_matrix(i,j) = params[k++];
235 }
236 }
237
238 template <unsigned int Tdim, unsigned int Fdim>
239 void
240 rgrl_est_proj_func<Tdim, Fdim>::
f(vnl_vector<double> const & x,vnl_vector<double> & fx)241 f(vnl_vector<double> const& x, vnl_vector<double>& fx)
242 {
243 typedef rgrl_match_set::const_from_iterator FIter;
244 typedef FIter::to_iterator TIter;
245
246 vnl_vector_fixed<double, Tdim> mapped;
247 vnl_matrix_fixed<double, Tdim,Tdim> error_proj_sqrt;
248 unsigned int ind = 0;
249
250 // retrieve the projection matrix
251 vnl_matrix_fixed<double, Tdim+1, Fdim+1> proj;
252 restored_centered_proj( proj, x );
253
254 for ( unsigned ms = 0; ms<matches_ptr_->size(); ++ms )
255 if ( (*matches_ptr_)[ms] != nullptr ) // if pointer is valid
256 {
257 rgrl_match_set const& one_set = *((*matches_ptr_)[ms]);
258 for ( FIter fi=one_set.from_begin(); fi!=one_set.from_end(); ++fi )
259 {
260 // map from point
261 vnl_vector_fixed<double, Fdim> from = fi.from_feature()->location();
262 map_loc( mapped, proj, from );
263
264 for ( TIter ti=fi.begin(); ti!=fi.end(); ++ti ) {
265 vnl_vector_fixed<double, Tdim> to = ti.to_feature()->location();
266 error_proj_sqrt = ti.to_feature()->error_projector_sqrt();
267 double const wgt = std::sqrt(ti.cumulative_weight());
268 vnl_vector_fixed<double, Tdim> diff = error_proj_sqrt * (mapped - to);
269
270 // fill in
271 for ( unsigned i=0; i<Tdim; ++i )
272 fx(ind++) = wgt * diff[i];
273 }
274 }
275 }
276
277 // check
278 assert( ind == get_number_of_residuals() );
279 }
280
281 template <unsigned int Tdim, unsigned int Fdim>
282 void
283 rgrl_est_proj_func<Tdim, Fdim>::
gradf(vnl_vector<double> const & x,vnl_matrix<double> & jacobian)284 gradf(vnl_vector<double> const& x, vnl_matrix<double>& jacobian)
285 {
286 typedef rgrl_match_set::const_from_iterator FIter;
287 typedef FIter::to_iterator TIter;
288
289 assert( jacobian.rows() == get_number_of_residuals() );
290 assert( jacobian.cols()+1 == proj_size_ );
291
292 vnl_matrix_fixed<double, Tdim, Tdim> error_proj_sqrt;
293 vnl_matrix_fixed<double, Tdim, proj_size_-1> base_jac, jac;
294
295 // retrieve the projection matrix
296 vnl_matrix_fixed<double, Tdim+1, Fdim+1> proj;
297 restored_centered_proj( proj, x );
298
299 unsigned int ind = 0;
300 for ( unsigned ms = 0; ms<matches_ptr_->size(); ++ms )
301 if ( (*matches_ptr_)[ms] ) { // if pointer is valid
302
303 rgrl_match_set const& one_set = *((*matches_ptr_)[ms]);
304 for ( FIter fi=one_set.from_begin(); fi!=one_set.from_end(); ++fi )
305 {
306 // map from point
307 vnl_vector_fixed<double, Fdim> from = fi.from_feature()->location();
308
309 // jacobian computation
310 reduced_proj_jacobian( base_jac, proj, from );
311
312 for ( TIter ti=fi.begin(); ti!=fi.end(); ++ti ) {
313 //vnl_double_2 to = ti.to_feature()->location();
314 error_proj_sqrt = ti.to_feature()->error_projector_sqrt();
315 double const wgt = std::sqrt(ti.cumulative_weight());
316 jac = error_proj_sqrt * base_jac;
317 jac *= wgt;
318
319 // fill in
320 for ( unsigned i=0; i<Tdim; ++i,++ind ) {
321 for ( unsigned j=0; j+1<proj_size_; ++j )
322 jacobian(ind,j) = jac(i, j);
323 }
324 }
325 }
326 }
327 }
328
329 template <unsigned int Tdim, unsigned int Fdim>
330 bool
331 rgrl_est_proj_func<Tdim, Fdim>::
projective_estimate(vnl_matrix_fixed<double,Tdim+1,Fdim+1> & proj,vnl_matrix<double> & full_covar,vnl_vector_fixed<double,Fdim> & from_centre,vnl_vector_fixed<double,Tdim> & to_centre)332 projective_estimate( vnl_matrix_fixed<double, Tdim+1, Fdim+1>& proj,
333 vnl_matrix<double>& full_covar,
334 vnl_vector_fixed<double, Fdim>& from_centre,
335 vnl_vector_fixed<double, Tdim>& to_centre )
336 {
337 // compute weighted centres
338 // this function is going to resize the vnl_vector, use temporary ones instead
339 vnl_vector<double> fc, tc;
340 if ( !rgrl_est_compute_weighted_centres( *matches_ptr_, fc, tc ) ) {
341 return false;
342 }
343 assert( fc.size() == from_centre.size() );
344 assert( tc.size() == to_centre.size() );
345 from_centre = fc;
346 to_centre = tc;
347
348 // convert parameters
349 vnl_vector<double> p;
350 this->convert_parameters( p, proj, from_centre, to_centre );
351
352 vnl_levenberg_marquardt lm( *this );
353 //lm.set_trace( true );
354 //lm.set_check_derivatives( 1 );
355 // we don't need it to be super accurate
356 lm.set_f_tolerance( relative_threshold_ );
357 lm.set_max_function_evals( max_num_iterations_ );
358
359 bool ret;
360 if ( has_gradient() )
361 ret = lm.minimize_using_gradient(p);
362 else
363 ret = lm.minimize_without_gradient(p);
364 if ( !ret ) {
365 std::cerr << "Levenberg-Marquardt failed\n";
366 lm.diagnose_outcome(std::cerr);
367 return false;
368 }
369 //lm.diagnose_outcome(std::cout);
370
371 // convert parameters back into matrix form
372 this->restored_centered_proj( proj, p );
373 //std::cout << "Final params=" << proj << std::endl;
374
375 // compute covariance
376 // Jac^\top * Jac is INVERSE of the covariance
377 //
378 const unsigned int proj_size_ = (Fdim+1)*(Tdim+1);
379 vnl_matrix<double> jac(get_number_of_residuals(), proj_size_-1);
380 this->gradf( p, jac );
381 //
382
383 // SVD decomposition:
384 // Jac = U W V^\top
385 vnl_svd<double> svd( jac, zero_svd_thres_ );
386 if ( svd.rank()+1 < proj_size_ ) {
387 std::cerr << "The covariance of projection matrix ranks less than "
388 << proj_size_-1 << "!\n"
389 << " The singular values are " << svd.W() << std::endl;
390 return false;
391 }
392
393 //invert the W matrix and square it
394 vnl_diag_matrix<double> invW( proj_size_-1 );
395 for ( unsigned i=0; i+1 < proj_size_; ++i )
396 invW[i] = vnl_math::sqr( 1.0/svd.W(i) );
397
398 //compute inverse of Jac^\top * Jac
399 const vnl_matrix<double> covar( svd.V() * invW * vnl_transpose( svd.V() ) );
400
401 // convert the covar to full dof+1 matrix
402 // fill in the row and column of the fixed element
403 // with 0
404 full_covar.set_size( proj_size_, proj_size_ );
405 full_covar.fill( 0.0 );
406
407 const unsigned int param_index = index_row_*(Fdim+1) + index_col_;
408 for ( unsigned i=0,ic=0; i<proj_size_; ++i,++ic ) {
409 if ( i==param_index ) { --ic; continue; }
410 for ( unsigned j=0,jc=0; j<proj_size_; ++j,++jc ) {
411 if ( j==param_index ) { --jc; continue; }
412 full_covar( i, j ) = covar( ic, jc );
413 }
414 }
415
416 return true;
417 }
418
419 #if 0 // commented out
420 // --------------------------------------------------------------------
421 template <unsigned int Tdim, unsigned int Fdim>
422 rgrl_est_proj_lm::
423 rgrl_est_proj_lm( bool with_grad )
424 : with_grad_( with_grad )
425 {
426 rgrl_estimator::set_param_dof( Fdim*Tdim-1 );
427
428 // default value
429 rgrl_nonlinear_estimator::set_max_num_iter( 50 );
430 rgrl_nonlinear_estimator::set_rel_thres( 1e-5 );
431 }
432
433 template <unsigned int Tdim, unsigned int Fdim>
434 bool
435 rgrl_est_proj_lm::
436 projective_estimate( vnl_matrix_fixed<double, Tdim+1, Fdim+1>& proj,
437 vnl_matrix<double>& full_covar,
438 vnl_vector_fixed<double, Fdim>& from_centre,
439 vnl_vector_fixed<double, Tdim>& to_centre,
440 rgrl_set_of<rgrl_match_set_sptr> const& matches ) const
441 {
442 // count the number of constraints/residuals
443 typedef rgrl_match_set::const_from_iterator FIter;
444 typedef FIter::to_iterator TIter;
445 unsigned int tot_num = 0;
446 for ( unsigned ms = 0; ms<matches.size(); ++ms )
447 if ( matches[ms] ) { // if pointer is valid
448
449 rgrl_match_set const& one_set = *(matches[ms]);
450 for ( FIter fi=one_set.from_begin(); fi!=one_set.from_end(); ++fi )
451 if ( fi.size() ) {
452 tot_num += fi.size() * fi.begin().to_feature()->dim(); // each point provides two constraints
453 }
454 }
455
456 // Determine the weighted centres for computing the more stable
457 // covariance matrix of homography parameters
458 //
459 if ( !compute_weighted_centres( matches, from_centre, to_centre ) )
460 return 0;
461 DebugMacro( 3, "From center: " << from_centre
462 <<" To center: " << to_centre << std::endl );
463
464 // construct least square cost function
465 rgrl_est_proj_func<Tdim, Fdim> proj_func( matches, tot_num, with_grad_ );
466
467 // convert parameters
468 vnl_vector<double> p;
469 proj_func.convert_parameters( p, proj, from_centre, to_centre );
470
471 vnl_levenberg_marquardt lm( proj_func );
472 //lm.set_trace( true );
473 //lm.set_check_derivatives( 2 );
474 // we don't need it to be super accurate
475 lm.set_f_tolerance( relative_threshold_ );
476 lm.set_max_function_evals( max_num_iterations_ );
477
478 bool ret;
479 if ( with_grad_ )
480 ret = lm.minimize_using_gradient(p);
481 else
482 ret = lm.minimize_without_gradient(p);
483 if ( !ret ) {
484 WarningMacro( "Levenberg-Marquardt failed" );
485 return false;
486 }
487 //lm.diagnose_outcome(std::cout);
488
489 // convert parameters back into matrix form
490 restored_centered_proj( proj, p );
491
492 // compute covariance
493 // Jac^\top * Jac is INVERSE of the covariance
494 //
495 const unsigned int proj_size_ = (Fdim+1)*(Tdim+1);
496 vnl_matrix<double> jac(tot_num, proj_size_-1);
497 proj_func.gradf( p, jac );
498 //
499
500 // SVD decomposition:
501 // Jac = U W V^\top
502 vnl_svd<double> svd( jac, 1e-4 );
503 if ( svd.rank()+1 < proj_size_ ) {
504 WarningMacro( "The covariance of homography ranks less than 8! ");
505 return 0;
506 }
507
508 //invert the W matrix and square it
509 vnl_diag_matrix<double> invW( proj_size_-1 );
510 for ( unsigned i=0; i+1<proj_size_; ++i )
511 invW[i] = vnl_math::sqr( 1.0/svd.W(i) );
512
513 //compute inverse of Jac^\top * Jac
514 const vnl_matrix<double> covar( svd.V() * invW * vnl_transpose( svd.V() ) );
515
516 // convert the covar to full dof+1 matrix
517 // fill in the row and column of the fixed element
518 // with 0
519 full_covar.set_size( proj_size_, proj_size_ );
520
521 const unsigned int param_index = index_row_*(Fdim+1) + index_col_;
522 for ( unsigned i=0,ic=0; i<proj_size_; ++i,++ic ) {
523 if ( i==param_index ) { --ic; continue; }
524 for ( unsigned j=0;jc=0; j<proj_size_; ++j,++jc ) {
525 if ( j==param_index ) { --jc; continue; }
526 full_covar( i, j ) = covar( ic, jc );
527 }
528 }
529 }
530 #endif // 0
531
532 #undef RGRL_EST_PROJ_FUNC_INSTANTIATE
533 #define RGRL_EST_PROJ_FUNC_INSTANTIATE( tdim, fdim ) \
534 template class rgrl_est_proj_func< tdim, fdim >
535
536 #endif //rgrl_est_proj_func_txx_
537