1 /** @file gsCurvatureSmoothing.h
2 
3     @brief Computes a closed B-spline curve with a smaller number of curvature
4     extrema compared to a given closed B-spline curve  i.e. some kind of
5     smoothing the curvature of the curve. This smoothting can be done with the
6     help of two methods - total variation and Hadenfeld's algorithm (see Jan
7     Hadenfeld, Iteratives Glätten von B-Spline Kurven und B-Spline Flächen,
8     Shaker Verlag, PhD Thesis)
9 
10     This file is part of the G+Smo library.
11 
12     This Source Code Form is subject to the terms of the Mozilla Public
13     License, v. 2.0. If a copy of the MPL was not distributed with this
14     file, You can obtain one at http://mozilla.org/MPL/2.0/.
15 
16     Author(s): M. Kapl
17 */
18 
19 #pragma once
20 
21 #include <iostream>
22 #include <stdexcept>
23 
24 #include <gsCore/gsGeometry.h>
25 #include <gsNurbs/gsBSpline.h>
26 #include <gsCore/gsLinearAlgebra.h>
27 
28 namespace gismo
29 {
30 /**
31     @brief Class for computing a closed B-spline curve with a smaller number of curvature extrema compared to a given closed B-spline curve
32 
33     i.e. some
34     kind of smoothing the curvature of the curve with the help of two different methods
35 
36     \ingroup Modeling
37 */
38 
39 template<class T>
40 class gsCurvatureSmoothing
41 {
42 public:
43     /// default constructor
gsCurvatureSmoothing()44     gsCurvatureSmoothing()
45     {
46         m_curve_original = NULL;
47         m_curve_smooth = NULL;
48     }
49 
50     /// constructor
gsCurvatureSmoothing(const gsBSpline<T> & init_curve,const gsMatrix<T> & param_values,const gsMatrix<T> & points)51     gsCurvatureSmoothing(const gsBSpline<T> & init_curve, const gsMatrix<T>& param_values, const gsMatrix<T>& points)
52     {
53         m_curve_original = &init_curve;
54         m_curve_smooth = init_curve.clone().release();
55         m_param_values = param_values;
56         m_points = points;
57     }
58 
59     /// Destructor
~gsCurvatureSmoothing()60     ~gsCurvatureSmoothing()
61     {
62         delete m_curve_smooth; // deleting the B-spline curve
63     }
64 
65 
66     /// gives back the original B-spline curve
curveOriginal()67     const gsBSpline<T> & curveOriginal() const { return *m_curve_original; }
68     /// gives back the smoother B-spline curve
curveSmooth()69     const gsBSpline<T> & curveSmooth() const { return *m_curve_smooth; }
70 
71 
72     /// smooth the curve by total variation -- computes the stepsize by itsself (with the help of a backtracking line search method)
73     /// this method should be used instead of the two methods void smoothTotalVariationSelectLamda
74     void smoothTotalVariation(const T omega1, const T omega2, const T lamda, const T tau, const unsigned iter=50);
75 
76     /// smooth the curve by total variation -- uses different stepsizes (in listlamdas) in the gradient descent method (look for the best from the list!)
77     /// if possible use the method void smoothTotalVariation
78     void smoothTotalVariationSelectLamda(const T omega1, const T omega2, const gsMatrix<T> listlamdas, const unsigned iter=50);
79 
80     /// smooth the curve by total variation -- uses always the same stepsize lamda - but which has to be chosen!!
81     /// if possible use the method void smoothTotalVariation
82     void smoothTotalVariationSelectLamda(const T omega1, const T omega2, const T lamda, const unsigned iter=50);
83 
84     /// smooth the curve by smoothing only one cofficient in each step using the Hadenfeld algorithm --- the usual Hadenfeld algorithm --this method should be used
85     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);
86 
87 
88     /// smooth the curve in one step for all coefficients using the Hadenfeld algorithm. Be aware of the fact
89     /// that it is not ensured that we get a nice result --- can work but do not have to work (not sure that it will converge!!)
90     /// if possible use method void smoothHadenfeld
91     void smoothAllHadenfeld(const unsigned smooth_degree=4, const unsigned iter=500);
92 
93     /// writes the smooth curve to a file, which can be visualized in Mathematica (Name of Mathematica File VisualizationCurvatureSmoothing)
94     void write(std::ostream &os);
95 
96     /// computes approximation error of the smoother curve to the original point cloud
97     void computeApproxError(T & error);
98 
99     /// computes the L^{2}-norm approximation error of the smoother curve
100     void computeApproxErrorL2(T & error);
101 
102     /// computes the L-max-norm approximation error of the smoother curve
103     void computeApproxErrorLMax(T & error);
104 
105     /// computes the max approximation error between the coeffeicientsof the original and the smoother curve
106     void computeApproxErrorCoef(T & error);
107 
108     /// computes the curvature error of the smoother curve
109     void computeCurvatureError(T & error);
110 
111 private:
112 
113     /// the original B-spline curve
114     const gsBSpline<T> * m_curve_original;
115     /// the smoother B-spline curve
116     gsBSpline<T> * m_curve_smooth;
117     /// the parameter values of the original point cloud
118     gsMatrix<T> m_param_values;
119     /// the points of the original point cloud
120     gsMatrix<T> m_points;
121 
122     /// computes all values and derivatives (up to three) at the parameter values u for the given coefs
123     void compute_AllValues(gsBSplineBasis<T> * basis, gsMatrix<T> u, gsMatrix<T> *coefs, gsMatrix<T> & values0, gsMatrix<T> & values1, gsMatrix<T> & values2, gsMatrix<T> & values3);
124 
125     /// computes the objective function for given coefs and omega1 and omega2 -- objective function = omega1*ApproximationFunction + omega2*CurvatureFunction
126     void compute_ObjectiveFunction(gsBSplineBasis<T> * basis, gsMatrix<T> *coefs, const T omega1, const T omega2, T &value);
127 
128     /// set the smooth curve to the the original curve
129     // TODO: What is exactly the purpose of this function; why do we want to change output?
reset(gsBSpline<T> * newCurve)130     void reset(gsBSpline<T> * newCurve)
131     {
132         if (m_curve_smooth)
133             delete m_curve_smooth;
134         m_curve_smooth = newCurve;
135     }
136 
137 }; // class gsCurvatureSmoothing
138 
139 } // namespace gismo
140 
141 #ifndef GISMO_BUILD_LIB
142 #include GISMO_HPP_HEADER(gsCurvatureSmoothing.hpp)
143 #endif
144