30 index_t m_degree=m_curve_smooth->degree();
36 index_t num_rows=current_coefs.rows()-m_degree;
37 index_t num_cols=current_coefs.cols();
68 T cond11,cond12,cond21,cond22,helpcond;
72 compute_ObjectiveFunction(basis,¤t_coefs,omega1,omega2,m_value0);
76 for(
unsigned i=0;i<iter;i++){
79 for(
index_t j=0;j<num_rows;j++){
80 for(
index_t k=0;k<num_cols;k++){
81 deriv_coefs1 = current_coefs;
82 deriv_coefs1(j,k)+=delta;
83 deriv_coefs2 = current_coefs;
84 deriv_coefs2(j,k)-=delta;
88 deriv_coefs1(j+num_rows,k)+=delta;
89 deriv_coefs2(j+num_rows,k)-=delta;
93 compute_ObjectiveFunction(basis,&deriv_coefs1,omega1,omega2,m_value1);
94 compute_ObjectiveFunction(basis,&deriv_coefs2,omega1,omega2,m_value2);
97 m_gradient(j,k)=(m_value1-m_value2)/((T)(2)*delta);
109 while( (cond11 > cond12) || (cond21 > cond22) ) {
111 different_coefs=current_coefs;
112 different_coefs.conservativeResize(num_rows,num_cols);
113 different_coefs=different_coefs-m_lamda*m_gradient;
114 different_coefs.conservativeResize(num_rows+m_degree,num_cols);
115 for(
index_t k=0;k<m_degree;k++){
116 different_coefs.row(num_rows+k)=different_coefs.row(k);
119 compute_ObjectiveFunction(basis,&different_coefs,omega1,omega2,lamda_value);
123 for(
index_t j=0;j<num_rows;j++){
124 for(
index_t k=0;k<num_cols;k++){
125 deriv_coefs3 = different_coefs;
126 deriv_coefs3(j,k)+=delta;
127 deriv_coefs4 = different_coefs;
128 deriv_coefs4(j,k)-=delta;
132 deriv_coefs3(j+num_rows,k)+=delta;
133 deriv_coefs4(j+num_rows,k)-=delta;
137 compute_ObjectiveFunction(basis,&deriv_coefs3,omega1,omega2,m_value1);
138 compute_ObjectiveFunction(basis,&deriv_coefs4,omega1,omega2,m_value2);
141 m_gradient2(j,k)=(m_value1-m_value2)/((T)(2)*delta);
150 for(
index_t j=0;j<m_gradient.rows();j++){
151 for(
index_t k=0;k<m_gradient.cols();k++){
153 helpcond+=(T)(2)*m_gradient(j,k)*m_gradient(j,k);
156 helpcond+=m_gradient(j,k)*m_gradient(j,k);
160 cond12=m_lamda*c1*helpcond+m_value0;
164 for(
index_t j=0;j<m_gradient.rows();j++){
165 for(
index_t k=0;k<m_gradient.cols();k++){
167 cond21+=(T)(2)*m_gradient(j,k)*m_gradient2(j,k);
170 cond21+=m_gradient(j,k)*m_gradient2(j,k);
185 current_coefs.conservativeResize(num_rows,num_cols);
186 current_coefs=current_coefs-m_lamda*m_gradient;
187 current_coefs.conservativeResize(num_rows+m_degree,num_cols);
188 for(
index_t k=0;k<m_degree;k++){
189 current_coefs.row(num_rows+k)=current_coefs.row(k);
193 compute_ObjectiveFunction(basis,¤t_coefs,omega1,omega2,m_value0);
195 gsDebug <<
"Step: " << i+1 <<
" lamda: " << m_lamda <<
" objective value: " <<m_value0 <<
"\n";
206 index_t m_degree=m_curve_smooth->degree();
212 index_t num_rows=current_coefs.rows()-m_degree;
213 index_t num_cols=current_coefs.cols();
237 compute_ObjectiveFunction(basis,¤t_coefs,omega1,omega2,m_value0);
239 for(
unsigned i=0;i<iter;i++){
242 for(
index_t j=0;j<num_rows;j++){
243 for(
index_t k=0;k<num_cols;k++){
244 deriv_coefs1 = current_coefs;
245 deriv_coefs1(j,k)+=delta;
246 deriv_coefs2 = current_coefs;
247 deriv_coefs2(j,k)-=delta;
251 deriv_coefs1(j+num_rows,k)+=delta;
252 deriv_coefs2(j+num_rows,k)-=delta;
256 compute_ObjectiveFunction(basis,&deriv_coefs1,omega1,omega2,m_value1);
257 compute_ObjectiveFunction(basis,&deriv_coefs2,omega1,omega2,m_value2);
260 m_gradient(j,k)=(m_value1-m_value2)/((T)(2)*delta);
270 for(
index_t j=0;j<listlamdas.cols();j++){
272 different_coefs=current_coefs;
274 different_coefs.conservativeResize(num_rows,num_cols);
275 different_coefs=different_coefs-listlamdas(0,j)*m_gradient;
276 different_coefs.conservativeResize(num_rows+m_degree,num_cols);
277 for(
index_t k=0;k<m_degree;k++){
278 different_coefs.row(num_rows+k)=different_coefs.row(k);
282 compute_ObjectiveFunction(basis,&different_coefs,omega1,omega2,lamda_value);
285 if(lamda_value<max_value){
286 max_value=lamda_value;
287 m_lamda=listlamdas(0,j);
292 current_coefs.conservativeResize(num_rows,num_cols);
293 current_coefs=current_coefs-m_lamda*m_gradient;
294 current_coefs.conservativeResize(num_rows+m_degree,num_cols);
295 for(
index_t k=0;k<m_degree;k++){
296 current_coefs.row(num_rows+k)=current_coefs.row(k);
300 compute_ObjectiveFunction(basis,¤t_coefs,omega1,omega2,m_value0);
302 gsDebug <<
"Step: " << i+1 <<
" lamda: " << m_lamda <<
" objective value: " <<m_value0 <<
"\n";
316 index_t m_degree=m_curve_smooth->degree();
320 index_t num_rows=current_coefs.rows()-m_degree;
321 index_t num_cols=current_coefs.cols();
340 compute_ObjectiveFunction(basis,¤t_coefs,omega1,omega2,m_value0);
342 for(
unsigned i=0;i<iter;i++){
345 for(
index_t j=0;j<num_rows;j++){
346 for(
index_t k=0;k<num_cols;k++){
347 deriv_coefs1 = current_coefs;
348 deriv_coefs1(j,k)+=delta;
349 deriv_coefs2 = current_coefs;
350 deriv_coefs2(j,k)-=delta;
354 deriv_coefs1(j+num_rows,k)+=delta;
355 deriv_coefs2(j+num_rows,k)-=delta;
359 compute_ObjectiveFunction(basis,&deriv_coefs1,omega1,omega2,m_value1);
360 compute_ObjectiveFunction(basis,&deriv_coefs2,omega1,omega2,m_value2);
363 m_gradient(j,k)=(m_value1-m_value2)/((T)(2)*delta);
368 current_coefs.conservativeResize(num_rows,num_cols);
369 current_coefs=current_coefs-lamda*m_gradient;
370 current_coefs.conservativeResize(num_rows+m_degree,num_cols);
371 for(
index_t k=0;k<m_degree;k++){
372 current_coefs.row(num_rows+k)=current_coefs.row(k);
376 compute_ObjectiveFunction(basis,¤t_coefs,omega1,omega2,m_value0);
378 gsDebug <<
"Step: " << i+1 <<
" lamda: " << lamda <<
" objective value: " <<m_value0 <<
"\n";
392 index_t m_degree=m_curve_smooth->degree();
394 index_t num_rows=m_coefs.rows()-m_degree;
395 m_coefs.conservativeResize(num_rows,2);
398 index_t m_iter_total= min(iter_total,iter_step*num_rows);
403 m_coefs_original=m_curve_original->coefs();
404 m_coefs_original.conservativeResize(num_rows,2);
408 m_coefs_original=m_coefs;
413 m_iterated.setZero();
415 T max_value,coef0,coef1,current0,current1,dist_value,m_smooth1,m_smooth2,m_smooth3,m_smooth4;
419 if(smooth_degree==2){
420 m_smooth1=(43.0/95.0);
421 m_smooth2=(16.0/95.0);
422 m_smooth3=-(11.0/95.0);
423 m_smooth4=-(1.0/190.0);
425 else if (smooth_degree==4){
427 m_smooth2=-(2.0/5.0);
428 m_smooth3=(4.0/35.0);
429 m_smooth4=-(1.0/70.0);
432 m_smooth1=(17.0/25.0);
433 m_smooth2=-(4.0/25.0);
434 m_smooth3=-(1.0/25.0);
435 m_smooth4=(1.0/50.0);
440 for(
index_t j=0;j<m_iter_total;j++){
446 for(
index_t i=0;i<num_rows;i++){
447 if(m_iterated(i)<iter_step){
449 current0=m_smooth1*m_coefs((num_rows+i-1)%num_rows,0)+m_smooth1*m_coefs((i+1)%num_rows,0)+m_smooth2*m_coefs((num_rows+i-2)%num_rows,0)
450 +m_smooth2*m_coefs((i+2)%num_rows,0)+m_smooth3*m_coefs((num_rows+i-3)%num_rows,0)+m_smooth3*m_coefs((i+3)%num_rows,0)
451 +m_smooth4*m_coefs((num_rows+i-4)%num_rows,0)+m_smooth4*m_coefs((i+4)%num_rows,0);
452 current1=m_smooth1*m_coefs((num_rows+i-1)%num_rows,1)+m_smooth1*m_coefs((i+1)%num_rows,1)+m_smooth2*m_coefs((num_rows+i-2)%num_rows,1)
453 +m_smooth2*m_coefs((i+2)%num_rows,1)+m_smooth3*m_coefs((num_rows+i-3)%num_rows,1)+m_smooth3*m_coefs((i+3)%num_rows,1)
454 +m_smooth4*m_coefs((num_rows+i-4)%num_rows,1)+m_smooth4*m_coefs((i+4)%num_rows,1);
457 dist_value=sqrt((current0-m_coefs(i,0))*(current0-m_coefs(i,0))+(current1-m_coefs(i,1))*(current1-m_coefs(i,1)));
459 if(dist_value>max_value){
463 max_value=dist_value;
471 max_value=sqrt((coef0-m_coefs_original(index,0))*(coef0-m_coefs_original(index,0))+(coef1-m_coefs_original(index,1))*(coef1-m_coefs_original(index,1)));
473 m_coefs(index,0)= m_coefs_original(index,0)+(delta/max_value)*(coef0-m_coefs_original(index,0));
474 m_coefs(index,1)= m_coefs_original(index,1)+(delta/max_value)*(coef1-m_coefs_original(index,1));
478 m_coefs(index,0)=coef0;
479 m_coefs(index,1)=coef1;
484 m_coefs.conservativeResize(num_rows+m_degree,2);
485 for(
index_t k=0;k<m_degree;k++){
486 m_coefs.row(num_rows+k)=m_coefs.row(k);
491 m_knots=m_curve_smooth->knots();
502 index_t m_degree=m_curve_smooth->degree();
504 index_t num_rows=m_coefs.rows()-m_degree;
505 m_coefs.conservativeResize(num_rows,2);
510 T m_smooth1, m_smooth2, m_smooth3, m_smooth4;
512 if(smooth_degree==2){
513 m_smooth1=(43.0/95.0);
514 m_smooth2=(16.0/95.0);
515 m_smooth3=-(11.0/95.0);
516 m_smooth4=-(1.0/190.0);
518 else if(smooth_degree==4){
520 m_smooth2=-(2.0/5.0);
521 m_smooth3=(4.0/35.0);
522 m_smooth4=-(1.0/70.0);
525 m_smooth1=(17.0/25.0);
526 m_smooth2=-(4.0/25.0);
527 m_smooth3=-(1.0/25.0);
528 m_smooth4=(1.0/50.0);
532 for (
index_t i=0;i<num_rows;i++){
533 m_A(i,(num_rows+i-1)%num_rows)=m_smooth1;
534 m_A(i,(i+1)%num_rows)=m_smooth1;
535 m_A(i,(num_rows+i-2)%num_rows)=m_smooth2;
536 m_A(i,(i+2)%num_rows)=m_smooth2;
537 m_A(i,(num_rows+i-3)%num_rows)=m_smooth3;
538 m_A(i,(i+3)%num_rows)=m_smooth3;
539 m_A(i,(num_rows+i-4)%num_rows)=m_smooth4;
540 m_A(i,(i+4)%num_rows)=m_smooth4;
544 for(
unsigned k=0;k<iter;k++){
550 m_coefs.conservativeResize(num_rows+m_degree,2);
551 for(
index_t k=0;k<m_degree;k++){
552 m_coefs.row(num_rows+k)=m_coefs.row(k);
557 m_knots=m_curve_smooth->knots();
568 index_t num_rows=m_coefs.rows();
570 for(
index_t k=0;k<num_rows-1;k++){
571 os <<
"{" << m_coefs(k,0) <<
"," << m_coefs(k,1) <<
"},";
573 os <<
"{" << m_coefs(num_rows-1,0) <<
"," << m_coefs(num_rows-1,1) <<
"}}\n";
581 m_curve_smooth->eval_into(m_param_values.transpose(),results);
582 results.transposeInPlace();
585 for(
index_t k=0;k<m_points.rows();k++){
586 error+=math::pow(m_points(k,0)-results(k,0),2)+math::pow(m_points(k,1)-results(k,1),2);
595 computeApproxError(error);
596 error= math::sqrt( error / static_cast<T>(m_points.rows()) );
604 m_curve_smooth->eval_into(m_param_values.transpose(),results);
605 results.transposeInPlace();
608 for(
index_t k=0;k<m_points.rows();k++)
610 error= math::max(error,(T)math::sqrt(math::pow(m_points(k,0)-results(k,0),2)+math::pow(m_points(k,1)-results(k,1),2)));
618 const gsMatrix<T> & coefs_original=m_curve_original->coefs();
619 const gsMatrix<T> & coefs_smooth=m_curve_smooth->coefs();
624 for(
index_t k=0;k<coefs_original.rows();k++)
626 error= math::max(error, (coefs_original.row(k) - coefs_smooth.row(k)).norm() );
637 gsMatrix<T> & current_coefs=m_curve_smooth->coefs();
645 compute_ObjectiveFunction(basis,¤t_coefs,0,1,error);
653 std::vector<gsMatrix<T> > m_results;
660 values0.resize(coefs->cols(),u.cols());
661 values1.resize(coefs->cols(),u.cols());
662 values2.resize(coefs->cols(),u.cols());
663 values3.resize(coefs->cols(),u.cols());
673 for(
index_t i=0;i<u.cols();i++)
676 values0.col(i)+=coefs->row(actives(k,i))*m_results[0](k,i);
677 values1.col(i)+=coefs->row(actives(k,i))*m_results[1](k,i);
678 values2.col(i)+=coefs->row(actives(k,i))*m_results[2](k,i);
679 values3.col(i)+=coefs->row(actives(k,i))*m_results[3](k,i);
698 compute_AllValues(basis,m_param_values.transpose(),coefs,m_values0,m_values1,m_values2,m_values3);
702 for(
index_t i=0;i<m_param_values.rows();i++){
703 objective1+=math::pow(m_values0(0,i)-m_points(i,0),2)+math::pow(m_values0(1,i)-m_points(i,1),2);
705 objective2+=
math::abs( 6.0*(m_values1(1,i)*m_values2(0,i) - m_values1(0,i)*m_values2(1,i))*(m_values1(0,i)*m_values2(0,i) + m_values1(1,i)*m_values2(1,i)) +
706 (T)(2)*( (math::pow(m_values1(0,i),2)+math::pow(m_values1(1,i),2)) * ((-1.0)*m_values1(1,i)*m_values3(0,i)+ m_values1(0,i)*m_values3(1,i)) ) )/
707 ((T)(2)*math::pow( math::pow(m_values1(0,i),2)+math::pow(m_values1(1,i),2) ,2.5) );
709 objective2=objective2/((T)(0.0+m_param_values.rows()));
712 value=omega1*objective1+omega2*objective2;
#define gsDebug
Definition: gsDebug.h:61
void smoothTotalVariationSelectLamda(const T omega1, const T omega2, const gsMatrix< T > listlamdas, const unsigned iter=50)
Definition: gsCurvatureSmoothing.hpp:203
void computeCurvatureError(T &error)
computes the curvature error of the smoother curve
Definition: gsCurvatureSmoothing.hpp:632
void computeApproxErrorCoef(T &error)
computes the max approximation error between the coeffeicientsof the original and the smoother curve ...
Definition: gsCurvatureSmoothing.hpp:616
S give(S &x)
Definition: gsMemory.h:266
void compute_ObjectiveFunction(gsBSplineBasis< T > *basis, gsMatrix< T > *coefs, const T omega1, const T omega2, T &value)
computes the objective function for given coefs and omega1 and omega2 – objective function = omega1*A...
Definition: gsCurvatureSmoothing.hpp:686
#define index_t
Definition: gsConfig.h:32
void computeApproxErrorL2(T &error)
computes the L^{2}-norm approximation error of the smoother curve
Definition: gsCurvatureSmoothing.hpp:591
A B-spline function of one argument, with arbitrary target dimension.
Definition: gsBSpline.h:50
void smoothTotalVariation(const T omega1, const T omega2, const T lamda, const T tau, const unsigned iter=50)
Definition: gsCurvatureSmoothing.hpp:27
Computes a closed B-spline curve with a smaller number of curvature extrema compared to a given close...
A univariate B-spline basis.
Definition: gsBSplineBasis.h:694
virtual void evalAllDers_into(const gsMatrix< T > &u, int n, std::vector< gsMatrix< T > > &result, bool sameElement=false) const
Evaluate the nonzero functions and their derivatives up to order n at points u into result...
Definition: gsBSplineBasis.hpp:864
void computeApproxErrorLMax(T &error)
computes the L-max-norm approximation error of the smoother curve
Definition: gsCurvatureSmoothing.hpp:601
void smoothHadenfeld(const unsigned smooth_degree, const T delta, const index_t iter_step, const index_t iter_total, gsVector< index_t > &iterated, const bool original=true)
smooth the curve by smoothing only one cofficient in each step using the Hadenfeld algorithm — the us...
Definition: gsCurvatureSmoothing.hpp:388
void compute_AllValues(gsBSplineBasis< T > *basis, gsMatrix< T > u, gsMatrix< T > *coefs, gsMatrix< T > &values0, gsMatrix< T > &values1, gsMatrix< T > &values2, gsMatrix< T > &values3)
computes all values and derivatives (up to three) at the parameter values u for the given coefs ...
Definition: gsCurvatureSmoothing.hpp:650
void smoothAllHadenfeld(const unsigned smooth_degree=4, const unsigned iter=500)
Definition: gsCurvatureSmoothing.hpp:500
void active_into(const gsMatrix< T > &u, gsMatrix< index_t > &result) const
Returns the indices of active basis functions at points u, as a list of indices, in result...
Definition: gsBSplineBasis.hpp:127
Class for representing a knot vector.
Definition: gsKnotVector.h:79
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
void write(std::ostream &os)
writes the smooth curve to a file, which can be visualized in Mathematica (Name of Mathematica File V...
Definition: gsCurvatureSmoothing.hpp:565
void computeApproxError(T &error)
computes approximation error of the smoother curve to the original point cloud
Definition: gsCurvatureSmoothing.hpp:578