diff --git a/Interpolation/include/CGAL/interpolation_functions.h b/Interpolation/include/CGAL/interpolation_functions.h index 535816258a8..879d5284cc4 100644 --- a/Interpolation/include/CGAL/interpolation_functions.h +++ b/Interpolation/include/CGAL/interpolation_functions.h @@ -32,20 +32,21 @@ #include namespace CGAL { - -//Functor class for accessing the function values/gradients + +// Functor class for accessing the function values/gradients template< class Map > struct Data_access - : public CGAL::unary_function > + : public CGAL::unary_function > { typedef typename Map::mapped_type Data_type; typedef typename Map::key_type Key_type; - Data_access< Map >(const Map& m): map(m){} + Data_access(const Map& m): map(m){} std::pair< Data_type, bool> - operator()(const Key_type& p) const { + operator()(const Key_type& p) const + { typename Map::const_iterator mit = map.find(p); if(mit!= map.end()) return std::make_pair(mit->second, true); @@ -56,7 +57,7 @@ struct Data_access }; //the interpolation functions: -template < class ForwardIterator, class Functor> +template < class ForwardIterator, class Functor > typename Functor::result_type::first_type linear_interpolation(ForwardIterator first, ForwardIterator beyond, const typename @@ -68,7 +69,8 @@ linear_interpolation(ForwardIterator first, ForwardIterator beyond, typedef typename Functor::result_type::first_type Value_type; Value_type result(0); typename Functor::result_type val; - for(; first !=beyond; ++first){ + for(; first !=beyond; ++first) + { val = function_value(first->first); CGAL_assertion(val.second); result += (first->second/norm) * val.first; @@ -77,7 +79,8 @@ linear_interpolation(ForwardIterator first, ForwardIterator beyond, } - template < class ForwardIterator, class Functor, class GradFunctor, class Traits, class Point> + template < class ForwardIterator, class Functor, class GradFunctor, + class Traits, class Point > std::pair< typename Functor::result_type::first_type, bool> quadratic_interpolation(ForwardIterator first, ForwardIterator beyond, const typename @@ -88,21 +91,22 @@ quadratic_interpolation(ForwardIterator first, ForwardIterator beyond, GradFunctor function_gradient, const Traits& traits) { - CGAL_precondition(norm >0); - Interpolation::internal::V2P v2p(traits); + CGAL_precondition(norm > 0); typedef typename Functor::result_type::first_type Value_type; + + Interpolation::internal::V2P v2p(traits); Value_type result(0); typename Functor::result_type f; typename GradFunctor::result_type grad; - for(; first !=beyond; ++first){ + for(; first !=beyond; ++first) + { f = function_value(first->first); grad = function_gradient(first->first); //test if value and gradient are correctly retrieved: CGAL_assertion(f.second); if(!grad.second) return std::make_pair(Value_type(0), false); - result += (first->second/norm) - *( f.first + grad.first* + result += (first->second/norm) * (f.first + grad.first * traits.construct_scaled_vector_d_object() (traits.construct_vector_d_object()(v2p(first->first), p),0.5)); } @@ -110,7 +114,8 @@ quadratic_interpolation(ForwardIterator first, ForwardIterator beyond, } - template < class ForwardIterator, class Functor, class GradFunctor, class Traits, class Point> + template < class ForwardIterator, class Functor, class GradFunctor, + class Traits, class Point > std::pair< typename Functor::result_type::first_type, bool> sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, const typename @@ -122,69 +127,68 @@ sibson_c1_interpolation(ForwardIterator first, ForwardIterator beyond, const Traits& traits) { CGAL_precondition(norm >0); - Interpolation::internal::V2P v2p(traits); typedef typename Functor::result_type::first_type Value_type; typedef typename Traits::FT Coord_type; Coord_type term1(0), term2(term1), term3(term1), term4(term1); - Value_type linear_int(0),gradient_int(0); + Value_type linear_int(0), gradient_int(0); typename Functor::result_type f; typename GradFunctor::result_type grad; + Interpolation::internal::V2P v2p(traits); - for(; first !=beyond; ++first){ + for(; first !=beyond; ++first) + { f = function_value(first->first); grad = function_gradient(first->first); CGAL_assertion(f.second); if(!grad.second) - //the values are not correct: - return std::make_pair(Value_type(0), false); + return std::make_pair(Value_type(0), false); //the values are not correct Coord_type coeff = first->second/norm; - Coord_type squared_dist = traits. - compute_squared_distance_d_object()(v2p(first->first), p); + Coord_type squared_dist = traits.compute_squared_distance_d_object()(v2p(first->first), p); Coord_type dist = CGAL_NTS sqrt(squared_dist); - if(squared_dist ==0){ + if(squared_dist == 0) + { ForwardIterator it = first; CGAL_USE(it); - CGAL_assertion(++it==beyond); + CGAL_assertion(++it == beyond); return std::make_pair(f.first, true); } + //three different terms to mix linear and gradient //interpolation - term1 += coeff / dist; - term2 += coeff * squared_dist; - term3 += coeff * dist; + term1 += coeff / dist; + term2 += coeff * squared_dist; + term3 += coeff * dist; linear_int += coeff * f.first; - gradient_int += (coeff/dist) - * (f.first + grad.first * + gradient_int += (coeff/dist) * (f.first + grad.first * traits.construct_vector_d_object()(v2p(first->first), p)); } - term4 = term3/ term1; + term4 = term3 / term1; gradient_int = gradient_int / term1; - return std::make_pair((term4* linear_int + term2 * gradient_int)/ + return std::make_pair((term4 * linear_int + term2 * gradient_int) / (term4 + term2), true); } -//this method works with rational number types: -//modification of Sibson's interpolant without sqrt -//following a proposition by Gunther Rote: +// This method works with rational number types: +// modification of Sibson's interpolant without sqrt +// following a proposition by Gunther Rote: // -// the general scheme: +// The general scheme: // Coord_type inv_weight = f(dist); //i.e. dist^2 // term1 += coeff/inv_weight; -// term2 += coeff * squared_dist; -// term3 += coeff*(squared_dist/inv_weight); -// gradient_int += (coeff/inv_weight)* -// (vh->get_value()+ vh->get_gradient() -// *(p - vh->point())); +// term2 += coeff * squared_dist; +// term3 += coeff*(squared_dist/inv_weight); +// gradient_int += (coeff/inv_weight) * (vh->get_value()+ vh->get_gradient() * (p - vh->point())); - template < class ForwardIterator, class Functor, class GradFunctor, class Traits, class Point> -std::pair< typename Functor::result_type::first_type, bool> +template < class ForwardIterator, class Functor, class GradFunctor, + class Traits, class Point > +std::pair< typename Functor::result_type::first_type, bool > sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond, const typename std::iterator_traits:: @@ -194,43 +198,46 @@ sibson_c1_interpolation_square(ForwardIterator first, ForwardIterator beyond, GradFunctor function_gradient, const Traits& traits) { - CGAL_precondition(norm >0); + CGAL_precondition(norm > 0); Interpolation::internal::V2P v2p(traits); typedef typename Functor::result_type::first_type Value_type; typedef typename Traits::FT Coord_type; Coord_type term1(0), term2(term1), term3(term1), term4(term1); - Value_type linear_int(0),gradient_int(0); + Value_type linear_int(0), gradient_int(0); typename Functor::result_type f; typename GradFunctor::result_type grad; - for(; first!=beyond; ++first){ + for(; first!=beyond; ++first) + { f = function_value(first->first); grad = function_gradient(first->first); CGAL_assertion(f.second); + if(!grad.second) - //the gradient is not known - return std::make_pair(Value_type(0), false); + return std::make_pair(Value_type(0), false); // the gradient is not known Coord_type coeff = first->second/norm; Coord_type squared_dist = traits.compute_squared_distance_d_object()(v2p(first->first), v2p(p)); - if(squared_dist ==0){ + if(squared_dist ==0) + { ForwardIterator it = first; CGAL_USE(it); - CGAL_assertion(++it==beyond); - return std::make_pair(f.first,true); + CGAL_assertion(++it == beyond); + return std::make_pair(f.first, true); } + //three different terms to mix linear and gradient //interpolation - term1 += coeff / squared_dist; - term2 += coeff * squared_dist; - term3 += coeff; + term1 += coeff / squared_dist; + term2 += coeff * squared_dist; + term3 += coeff; linear_int += coeff * f.first; gradient_int += (coeff/squared_dist) * (f.first + grad.first * - traits.construct_vector_d_object()(v2p(first->first), v2p(p))); + traits.construct_vector_d_object()(v2p(first->first), v2p(p))); } term4 = term3/ term1; @@ -264,9 +271,10 @@ farin_c1_interpolation(RandomAccessIterator first, typename Functor::result_type f; typename GradFunctor::result_type grad; - int n= static_cast(beyond - first); - if( n==1){ - f= function_value(first->first); + int n = static_cast(beyond - first); + if(n == 1) + { + f = function_value(first->first); CGAL_assertion(f.second); return std::make_pair(f.first, true); } @@ -279,31 +287,32 @@ farin_c1_interpolation(RandomAccessIterator first, Value_type result(0); const Coord_type fac3(3); - std::vector< std::vector > - ordinates(n,std::vector(n, Value_type(0))); + std::vector< std::vector > ordinates(n, std::vector(n, Value_type(0))); - for(int i =0; isecond); - //for later: the function value of it->first: + // for later: the function value of it->first: f = function_value(it->first); CGAL_assertion(f.second); ordinates[i][i] = f.first; - //control point = data point + // control point = data point result += coord_i_square * it->second* ordinates[i][i]; - //compute tangent plane control point (one 2, one 1 entry) + // compute tangent plane control point (one 2, one 1 entry) Value_type res_i(0); - for(int j =0; jfirst); if(!grad.second) - //the gradient is not known - return std::make_pair(Value_type(0), false); + return std::make_pair(Value_type(0), false); // the gradient is not known //ordinates[i][j] = (p_j - p_i) * g_i ordinates[i][j] = grad.first * @@ -315,33 +324,38 @@ farin_c1_interpolation(RandomAccessIterator first, res_i += (fac3 * ordinates[i][i] + ordinates[i][j])* it2->second; } } - //res_i already multiplied by three + // res_i already multiplied by three result += coord_i_square *res_i; } - //the third type of control points: three 1 entries i,j,k + // the third type of control points: three 1 entries i,j,k for(int i=0; i< n; ++i) + { for(int j=i+1; j< n; ++j) - for(int k=j+1; k 6 * b_ijk = 3*(f_i + f_j + f_k) + 0.5*a - result += (Coord_type(2.0)*( ordinates[i][i]+ ordinates[j][j]+ - ordinates[k][k]) + result += (Coord_type(2.0)*(ordinates[i][i]+ ordinates[j][j]+ + ordinates[k][k]) + Coord_type(0.5)*(ordinates[i][j] + ordinates[i][k] + ordinates[j][i] + ordinates[j][k] + ordinates[k][i]+ ordinates[k][j])) - *(first+i)->second *(first+j)->second *(first+k)->second ; + * (first+i)->second *(first+j)->second *(first+k)->second ; } + } + } return std::make_pair(result/(CGAL_NTS square(norm)*norm), true); } -} //namespace CGAL +} // namespace CGAL #endif // CGAL_INTERPOLATION_FUNCTIONS_H diff --git a/Interpolation/include/CGAL/sibson_gradient_fitting.h b/Interpolation/include/CGAL/sibson_gradient_fitting.h index 536b1853ce6..1cfa86cab45 100644 --- a/Interpolation/include/CGAL/sibson_gradient_fitting.h +++ b/Interpolation/include/CGAL/sibson_gradient_fitting.h @@ -24,73 +24,69 @@ #include -#include #include #include #include -#include -#include +#include #include #include +#include +#include +#include + namespace CGAL { -template < class ForwardIterator, class Functor, class Traits> +template < class ForwardIterator, class Functor, class Traits > typename Traits::Vector_d sibson_gradient_fitting(ForwardIterator first, ForwardIterator beyond, const typename - std::iterator_traits:: - value_type::second_type& norm, + std::iterator_traits::value_type::second_type& norm, const typename Traits::Point_d& bare_p, const typename Functor::result_type::first_type fn, Functor function_value, const Traits& traits) { + CGAL_precondition( first != beyond && norm != 0); + Interpolation::internal::V2P v2p(traits); - CGAL_precondition( first!=beyond && norm!=0); typedef typename Traits::Aff_transformation_d Aff_transformation; typedef typename Traits::FT Coord_type; - typename Traits::Vector_d pn = - traits.construct_vector_d_object()(NULL_VECTOR); - Aff_transformation scaling, m, - Hn(traits.construct_null_matrix_d_object()()); + typename Traits::Vector_d pn = traits.construct_vector_d_object()(NULL_VECTOR); + Aff_transformation scaling, m, Hn(traits.construct_null_matrix_d_object()()); - for(;first!=beyond; ++first) + for(; first!=beyond; ++first) { const typename Traits::Point_d& bare_f = v2p(first->first); Coord_type square_dist = traits.compute_squared_distance_d_object()(bare_f, bare_p); CGAL_assertion(square_dist != 0); - Coord_type scale = first->second / (norm*square_dist); + Coord_type scale = first->second / (norm * square_dist); typename Traits::Vector_d d = traits.construct_vector_d_object()(bare_p, bare_f); - //compute the vector pn: + // compute the vector pn: typename Functor::result_type f = function_value(first->first); - CGAL_assertion(f.second);//function value of first->first is valid - pn = pn + traits.construct_scaled_vector_d_object() - (d,scale * (f.first - fn)); + CGAL_assertion(f.second); //function value of first->first is valid + pn = pn + traits.construct_scaled_vector_d_object()(d, scale * (f.first - fn)); //compute the matrix Hn: m = traits.construct_outer_product_d_object()(d); scaling = traits.construct_scaling_matrix_d_object()(scale); - Hn = traits.construct_sum_matrix_d_object()(Hn, scaling * m); } return Hn.inverse().transform(pn); } -template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH> +template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH > typename Traits::Vector_d sibson_gradient_fitting(const Triangul& tr, ForwardIterator first, ForwardIterator beyond, - const typename - std::iterator_traits:: - value_type::second_type& norm, + const typename std::iterator_traits::value_type::second_type& norm, VH vh, Functor function_value, const Traits& traits, @@ -110,14 +106,13 @@ sibson_gradient_fitting(const Triangul& tr, traits); } -template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH> +template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH > typename Traits::Vector_d sibson_gradient_fitting(const Triangul& tr, ForwardIterator first, ForwardIterator beyond, const typename - std::iterator_traits:: - value_type::second_type& norm, + std::iterator_traits::value_type::second_type& norm, VH vh, Functor function_value, const Traits& traits, @@ -136,14 +131,13 @@ sibson_gradient_fitting(const Triangul& tr, traits); } -template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH> +template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH > typename Traits::Vector_d sibson_gradient_fitting(const Triangul& tr, ForwardIterator first, ForwardIterator beyond, const typename - std::iterator_traits:: - value_type::second_type& norm, + std::iterator_traits::value_type::second_type& norm, VH vh, Functor function_value, const Traits& traits, @@ -151,7 +145,7 @@ sibson_gradient_fitting(const Triangul& tr, { const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point()); typename Functor::result_type fn = function_value(vh); - CGAL_assertion(fn.second); + CGAL_assertion(fn.second); return sibson_gradient_fitting(first, beyond, @@ -164,14 +158,14 @@ sibson_gradient_fitting(const Triangul& tr, template < class Triangul, class OutputIterator, class Functor, - class CoordFunctor, class OIF, class Traits> + class CoordFunctor, class OIF, class Traits > OutputIterator sibson_gradient_fitting_internal(const Triangul& tr, - OutputIterator out, - Functor function_value, - CoordFunctor compute_coordinates, - OIF fct, - const Traits& traits) + OutputIterator out, + Functor function_value, + CoordFunctor compute_coordinates, + OIF fct, + const Traits& traits) { typedef typename Traits::FT Coord_type; @@ -205,12 +199,12 @@ sibson_gradient_fitting_internal(const Triangul& tr, } -//the following functions allow to fit the gradients for all points in +// The following functions allow to fit the gradients for all points in // a triangulation except the convex hull points. // -> _nn2: natural_neighbor_coordinates_2 // -> _rn2: regular_neighbor_coordinates_2 // -> _sn2_3: surface_neighbor_coordinates_2_3 -template < class Dt, class OutputIterator, class Functor, class OIF, class Traits> +template < class Dt, class OutputIterator, class Functor, class OIF, class Traits > OutputIterator sibson_gradient_fitting_nn_2(const Dt& dt, OutputIterator out, @@ -238,7 +232,7 @@ sibson_gradient_fitting_nn_2(const Dt& dt, } -template < class Dt, class OutputIterator, class Functor, class Traits> +template < class Dt, class OutputIterator, class Functor, class Traits > OutputIterator sibson_gradient_fitting_nn_2(const Dt& dt, OutputIterator out, @@ -250,7 +244,7 @@ sibson_gradient_fitting_nn_2(const Dt& dt, } -template < class Rt, class OutputIterator, class Functor, class OIF, class Traits> +template < class Rt, class OutputIterator, class Functor, class OIF, class Traits > OutputIterator sibson_gradient_fitting_rn_2(const Rt& rt, OutputIterator out, @@ -263,7 +257,6 @@ sibson_gradient_fitting_rn_2(const Rt& rt, std::pair< typename Functor::argument_type, typename Traits::FT > > > CoordInserter; - typedef typename boost::mpl::if_< boost::is_same, Interpolation::internal::Vertex2WPoint, @@ -271,14 +264,14 @@ sibson_gradient_fitting_rn_2(const Rt& rt, >::type Coord_OIF; return sibson_gradient_fitting_internal(rt, out, function_value, - regular_neighbor_coordinates_2_object(), - fct, - traits); + regular_neighbor_coordinates_2_object(), + fct, + traits); } -template < class Rt, class OutputIterator, class Functor, class Traits> +template < class Rt, class OutputIterator, class Functor, class Traits > OutputIterator sibson_gradient_fitting_rn_2(const Rt& rt, OutputIterator out,