Code-Eli  0.3.6
curvature.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef eli_geom_surface_curvature_hpp
14 #define eli_geom_surface_curvature_hpp
15 
16 #include <iomanip>
17 #include <limits> // numeric_limits
18 
19 #include "eli/code_eli.hpp"
20 
21 namespace eli
22 {
23  namespace geom
24  {
25  namespace surface
26  {
27  namespace internal
28  {
29  template<typename surface__>
30  void calculate_surface_terms(typename surface__::data_type &e, typename surface__::data_type &f, typename surface__::data_type &g,
31  typename surface__::data_type &ee, typename surface__::data_type &ff, typename surface__::data_type &gg,
32  const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
33  {
34  typename surface__::point_type Su, Sv, Suu, Suv, Svv, n;
35 
36  // surfaces in two dimensions have no curvature
37  if (Su.innerSize()==2)
38  {
39  e=0;
40  f=0;
41  g=0;
42  ee=1;
43  ff=1;
44  gg=1;
45  return;
46  }
47  else if (Su.innerSize()!=3)
48  {
49  // how generalize curvature to surfaces in hyper-space?
50  assert(false);
51  return;
52  }
53 
54  // calculate the surface derivatives
55  Su=s.f_u(u, v);
56  Sv=s.f_v(u, v);
57  Suu=s.f_uu(u, v);
58  Suv=s.f_uv(u, v);
59  Svv=s.f_vv(u, v);
60  n=s.normal(u, v);
61 
62  // calculate terms
63  e=n.dot(Suu);
64  f=n.dot(Suv);
65  g=n.dot(Svv);
66  ee=Su.dot(Su);
67  ff=Su.dot(Sv);
68  gg=Sv.dot(Sv);
69  assert((ee*gg-ff*ff)>=0);
70  assert(std::abs(1-std::sqrt(ee*gg-ff*ff)/Su.cross(Sv).norm())<1e-4);
71  }
72 
73  template<typename surface__>
74  void principal_curvature_calc(typename surface__::data_type &kmax, typename surface__::data_type &kmin,
75  typename surface__::point_type &kmax_dir, typename surface__::point_type &kmin_dir, typename surface__::point_type &n,
76  const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v, bool calc_dir)
77  {
78  typename surface__::data_type H, K, tmp;
79 
80  // check to make sure have valid curve
81  assert(s.degree_u()>0);
82  assert(s.degree_v()>0);
83 
84  // check to make sure given valid parametric value
85  assert((u>=0) && (u<=1));
86  assert((v>=0) && (v<=1));
87 
88  // calculate the surface parameters
89  typename surface__::data_type e, f, g, ee, ff, gg;
90 
91  internal::calculate_surface_terms(e, f, g, ee, ff, gg, s, u, v);
92 
93  // calculate the mean and gaussian curvatures
94  tmp=ee*gg-ff*ff;
95  if (tmp==0)
96  {
97  kmax=std::numeric_limits<typename surface__::data_type>::max();
98  kmin=std::numeric_limits<typename surface__::data_type>::max();
99  return;
100  }
101  else
102  {
103  H=(e*gg-2*f*ff+g*ee)/(2*tmp);
104  K=(e*g-f*f)/tmp;
105  }
106 
107 
108  // calculate the principal curvatures
109  tmp=std::sqrt(H*H-K);
110  kmax=H+tmp;
111  kmin=H-tmp;
112 
113  // check if round-off error has crept in
114  typename surface__::tolerance_type tol;
115  if ((tmp<1e-5) && (tmp/std::max(e*e*gg*gg, std::max(f*f*ff*ff, g*g*ee*ee))<tol.get_relative_tolerance()))
116  {
117  kmax=H;
118  kmin=H;
119  tmp=0;
120  }
121 
122  // calculate the principal curvature directions
123  if (calc_dir)
124  {
125  typename surface__::point_type Su, Sv;
126  typename surface__::data_type lambda_max, lambda_min;
127 
128  Su=s.f_u(u, v);
129  Sv=s.f_v(u, v);
130  n=Su.cross(Sv);
131  n.normalize();
132 
133  // handle umbilic point case specially
134  if (tmp==0)
135  {
136  kmax_dir=Su;
137  kmin_dir=n.cross(kmax_dir);
138  }
139  else
140  {
141  tmp=g-kmax*gg;
142  if (tmp==0)
143  {
144  kmax_dir=Sv;
145  }
146  else
147  {
148  lambda_max=-(f-kmax*ff)/tmp;
149  kmax_dir=Su+lambda_max*Sv;
150  }
151 
152  tmp=g-kmin*gg;
153  if (tmp==0)
154  {
155  kmin_dir=Sv;
156  }
157  else
158  {
159  lambda_min=-(f-kmin*ff)/tmp;
160  kmin_dir=Su+lambda_min*Sv;
161  }
162  }
163  kmax_dir.normalize();
164  kmin_dir.normalize();
165 
166  // need to return right-handed system kmax_dir:kmin_dir:n so check direction of kmin_dir
167  tmp=n.cross(kmax_dir).dot(kmin_dir);
168  if (tmp<0)
169  {
170  kmin_dir*=-1;
171  }
172  assert(std::abs(1-std::abs(tmp))<0.001); // tmp1 should be very close to 1 or -1
173  }
174  }
175  }
176 
177  template<typename surface__>
178  void mean_curvature(typename surface__::data_type &k, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
179  {
180  // check to make sure have valid curve
181  assert(s.degree_u()>0);
182  assert(s.degree_v()>0);
183 
184  // check to make sure given valid parametric value
185  assert((u>=0) && (u<=1));
186  assert((v>=0) && (v<=1));
187 
188  // calculate the surface parameters
189  typename surface__::data_type e, f, g, ee, ff, gg;
190 
191  internal::calculate_surface_terms(e, f, g, ee, ff, gg, s, u, v);
192 
193  // calculate the curvature
194  typename surface__::data_type tmp;
195  tmp=ee*gg-ff*ff;
196  if ((tmp)==0)
197  {
198  k=std::numeric_limits<typename surface__::data_type>::max();
199  return;
200  }
201 
202  k=(e*gg-2*f*ff+g*ee)/(2*tmp);
203  }
204 
205  template<typename surface__>
206  void gaussian_curvature(typename surface__::data_type &k, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
207  {
208  // check to make sure have valid curve
209  assert(s.degree_u()>0);
210  assert(s.degree_v()>0);
211 
212  // check to make sure given valid parametric value
213  assert((u>=0) && (u<=1));
214  assert((v>=0) && (v<=1));
215 
216  // calculate the surface parameters
217  typename surface__::data_type e, f, g, ee, ff, gg;
218 
219  internal::calculate_surface_terms(e, f, g, ee, ff, gg, s, u, v);
220 
221  // calculate the curvature
222  typename surface__::data_type tmp;
223  tmp=ee*gg-ff*ff;
224  if (tmp==0)
225  {
226  k=std::numeric_limits<typename surface__::data_type>::max();
227  return;
228  }
229 
230  k=(e*g-f*f)/tmp;
231  }
232 
233  template<typename surface__>
234  void principal_curvature(typename surface__::data_type &kmax, typename surface__::data_type &kmin, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
235  {
236  typename surface__::point_type kmax_dir, kmin_dir, n;
237  internal::principal_curvature_calc(kmax, kmin, kmax_dir, kmin_dir, n, s, u, v, false);
238  }
239 
240  template<typename surface__>
241  void principal_curvature(typename surface__::data_type &kmax, typename surface__::data_type &kmin,
242  typename surface__::point_type &kmax_dir, typename surface__::point_type &kmin_dir, typename surface__::point_type &n,
243  const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
244  {
245  internal::principal_curvature_calc(kmax, kmin, kmax_dir, kmin_dir, n, s, u, v, true);
246  }
247  }
248  }
249 }
250 
251 #endif
Definition: math.hpp:20
void gaussian_curvature(typename surface__::data_type &k, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
Definition: curvature.hpp:206
void principal_curvature_calc(typename surface__::data_type &kmax, typename surface__::data_type &kmin, typename surface__::point_type &kmax_dir, typename surface__::point_type &kmin_dir, typename surface__::point_type &n, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v, bool calc_dir)
Definition: curvature.hpp:74
void principal_curvature(typename surface__::data_type &kmax, typename surface__::data_type &kmin, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
Definition: curvature.hpp:234
void mean_curvature(typename surface__::data_type &k, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
Definition: curvature.hpp:178
void calculate_surface_terms(typename surface__::data_type &e, typename surface__::data_type &f, typename surface__::data_type &g, typename surface__::data_type &ee, typename surface__::data_type &ff, typename surface__::data_type &gg, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
Definition: curvature.hpp:30