BRL-CAD
nurb_c2.c
Go to the documentation of this file.
1 /* N U R B _ C 2 . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1990-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @addtogroup nurb */
21 /** @{ */
22 /** @file primitives/bspline/nurb_c2.c
23  *
24  * Given parametric u, v values, return the curvature of the surface.
25  *
26  */
27 /** @} */
28 
29 #include "common.h"
30 
31 #include <math.h>
32 #include "bio.h"
33 
34 #include "vmath.h"
35 #include "raytrace.h"
36 #include "nurb.h"
37 
38 
39 void
40 rt_nurb_curvature(struct curvature *cvp, const struct face_g_snurb *srf, fastf_t u, fastf_t v)
41 {
42  struct face_g_snurb * us, *vs, * uus, * vvs, *uvs;
43  fastf_t ue[4], ve[4], uue[4], vve[4], uve[4], se[4];
44  fastf_t E, F, G; /* First Fundamental Form */
45  fastf_t L, M, N; /* Second Fundamental form */
46  fastf_t denom;
47  fastf_t wein[4]; /*Weingarten matrix */
48  fastf_t evec[3];
49  fastf_t mean, gauss, discrim;
50  vect_t norm;
51  int i;
52 
53  us = rt_nurb_s_diff(srf, RT_NURB_SPLIT_ROW);
54  vs = rt_nurb_s_diff(srf, RT_NURB_SPLIT_COL);
55  uus = rt_nurb_s_diff(us, RT_NURB_SPLIT_ROW);
56  vvs = rt_nurb_s_diff(vs, RT_NURB_SPLIT_COL);
57  uvs = rt_nurb_s_diff(vs, RT_NURB_SPLIT_ROW);
58 
59  rt_nurb_s_eval(srf, u, v, se);
60  rt_nurb_s_eval(us, u, v, ue);
61  rt_nurb_s_eval(vs, u, v, ve);
62  rt_nurb_s_eval(uus, u, v, uue);
63  rt_nurb_s_eval(vvs, u, v, vve);
64  rt_nurb_s_eval(uvs, u, v, uve);
65 
66  rt_nurb_free_snurb(us, (struct resource *)NULL);
67  rt_nurb_free_snurb(vs, (struct resource *)NULL);
68  rt_nurb_free_snurb(uus, (struct resource *)NULL);
69  rt_nurb_free_snurb(vvs, (struct resource *)NULL);
70  rt_nurb_free_snurb(uvs, (struct resource *)NULL);
71 
72  if (RT_NURB_IS_PT_RATIONAL(srf->pt_type)) {
73  for (i = 0; i < 3; i++) {
74  ue[i] = (1.0 / se[3] * ue[i]) -
75  (ue[3]/se[3]) * se[0]/se[3];
76  ve[i] = (1.0 / se[3] * ve[i]) -
77  (ve[3]/se[3]) * se[0]/se[3];
78  }
79  VCROSS(norm, ue, ve);
80  VUNITIZE(norm);
81  E = VDOT(ue, ue);
82  F = VDOT(ue, ve);
83  G = VDOT(ve, ve);
84 
85  for (i = 0; i < 3; i++) {
86  uue[i] = (1.0 / se[3] * uue[i]) -
87  2 * (uue[3]/se[3]) * uue[i] -
88  uue[3]/se[3] * (se[i]/se[3]);
89 
90  vve[i] = (1.0 / se[3] * vve[i]) -
91  2 * (vve[3]/se[3]) * vve[i] -
92  vve[3]/se[3] * (se[i]/se[3]);
93 
94  uve[i] = 1.0 / se[3] * uve[i] +
95  (-1.0 / (se[3] * se[3])) *
96  (ve[3] * ue[i] + ue[3] * ve[i] +
97  uve[3] * se[i]) +
98  (-2.0 / (se[3] * se[3] * se[3])) *
99  (ve[3] * ue[3] * se[i]);
100  }
101 
102  L = VDOT(norm, uue);
103  M = VDOT(norm, uve);
104  N = VDOT(norm, vve);
105 
106  } else {
107  VCROSS(norm, ue, ve);
108  VUNITIZE(norm);
109  E = VDOT(ue, ue);
110  F = VDOT(ue, ve);
111  G = VDOT(ve, ve);
112 
113  L = VDOT(norm, uue);
114  M = VDOT(norm, uve);
115  N = VDOT(norm, vve);
116  }
117 
118  if (srf->order[0] <= 2 && srf->order[1] <= 2) {
119  cvp->crv_c1 = cvp->crv_c2 = 0;
120  bn_vec_ortho(cvp->crv_pdir, norm);
121  return;
122  }
123 
124  denom = ((E*G) - (F*F));
125  gauss = (L * N - M *M)/denom;
126  mean = (G * L + E * N - 2 * F * M) / (2 * denom);
127  discrim = sqrt(mean * mean - gauss);
128 
129  cvp->crv_c1 = mean - discrim;
130  cvp->crv_c2 = mean + discrim;
131 
132  if (fabs(E*G - F*F) < 0.0001) {
133  /* XXX */
134  bu_log("rt_nurb_curvature: first fundamental form is singular E = %g F= %g G = %g\n",
135  E, F, G);
136  bn_vec_ortho(cvp->crv_pdir, norm); /* sanity */
137  return;
138  }
139 
140  wein[0] = ((G * L) - (F * M))/ (denom);
141  wein[1] = ((G * M) - (F * N))/ (denom);
142  wein[2] = ((E * M) - (F * L))/ (denom);
143  wein[3] = ((E * N) - (F * M))/ (denom);
144 
145  if (fabs(wein[1]) < 0.0001 && fabs(wein[3] - cvp->crv_c1) < 0.0001) {
146  evec[0] = 0.0; evec[1] = 1.0;
147  } else {
148  evec[0] = 1.0;
149  if (fabs(wein[1]) > fabs(wein[3] - cvp->crv_c1)) {
150  evec[1] = (cvp->crv_c1 - wein[0]) / wein[1];
151  } else {
152  evec[1] = wein[2] / (cvp->crv_c1 - wein[3]);
153  }
154  }
155 
156  cvp->crv_pdir[0] = evec[0] * ue[0] + evec[1] * ve[0];
157  cvp->crv_pdir[1] = evec[0] * ue[1] + evec[1] * ve[1];
158  cvp->crv_pdir[2] = evec[0] * ue[2] + evec[1] * ve[2];
159  VUNITIZE(cvp->crv_pdir);
160 }
161 
162 
163 /*
164  * Local Variables:
165  * mode: C
166  * tab-width: 8
167  * indent-tabs-mode: t
168  * c-file-style: "stroustrup"
169  * End:
170  * ex: shiftwidth=4 tabstop=8
171  */
void rt_nurb_free_snurb(struct face_g_snurb *srf, struct resource *res)
Definition: nurb_util.c:127
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
vect_t crv_pdir
Principle direction.
Definition: raytrace.h:307
#define N
Definition: randmt.c:39
void rt_nurb_curvature(struct curvature *cvp, const struct face_g_snurb *srf, fastf_t u, fastf_t v)
Definition: nurb_c2.c:40
Header file for the BRL-CAD common definitions.
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
void bn_vec_ortho(vect_t out, const vect_t in)
void rt_nurb_s_eval(const struct face_g_snurb *srf, fastf_t u, fastf_t v, fastf_t *final_value)
Definition: nurb_eval.c:49
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
struct face_g_snurb * rt_nurb_s_diff(const struct face_g_snurb *srf, int dir)
Definition: nurb_diff.c:58
double fastf_t
Definition: defines.h:300
#define M
Definition: msr.c:52