BRL-CAD
dspline.c
Go to the documentation of this file.
1 /* D S P L I N E . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2004-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 librt */
21 /** @{ */
22 /** @file librt/dspline.c
23  *
24  * Simple data (double) spline package.
25  *
26  * rt_dspline_matrix(m, type, tension, bias) create basis matrix
27  * rt_dspline4(m, a, b, c, d, alpha) interpolate 1 value
28  * rt_dspline4v(m, a, b, c, d, depth alpha) interpolate vectors
29  * rt_dspline(r, m, knots, n, depth, alpha) interpolate n knots over 0..1
30  *
31  * Example:
32  *
33  * mat_t m;
34  * double d;
35  * vect_t v;
36  * vect_t kn = { {0., 0., 0.},
37  * {0., 1., 0.},
38  * {.5, 1., 0.},
39  * {1., 1., 0.},
40  * {1., 0., 0.} };
41  *
42  * rt_dspline_matrix(m, "Beta", 0.5, 1.0);
43  *
44  * d = rt_dspline4(m, .0, .0, 1.0, 1.0, 0.25);
45  *
46  * for (p = 0.0; p <= 1.0; p += 0.0625) {
47  * rt_dspline(v, m, kn, 5, 3, p);
48  * bu_log("%g (%g %g %g)\n", p, V3ARGS(v));
49  * }
50  *
51  */
52 /** @} */
53 
54 #include "common.h"
55 
56 #include <stdlib.h>
57 #include <string.h>
58 #include <math.h>
59 #include "bio.h"
60 
61 #include "vmath.h"
62 
63 #include "raytrace.h"
64 
65 
66 static void
67 GetBeta(fastf_t *m, const double bias, const double tension)
68 {
69  register int i;
70  double d, b2, b3;
71  register double tmp;
72 
73  b2 = bias * bias;
74  b3 = bias * b2;
75 
76  tmp = 2.0 * b3;
77  m[12] = tmp;
78  m[ 0] = -tmp;
79 
80  tmp = tension+b2+bias;
81  m[ 1] = 2.0 * (tmp + b3);
82  m[ 2] = -2.0 * (tmp + 1.0);
83 
84  tmp = tension + 2.0 * b2;
85  m[ 5] = -3.0 * (tmp + 2.0 * b3);
86  m[ 6] = 3.0 * tmp;
87 
88  tmp = 6.0 * b3;
89  m[ 4] = tmp;
90  m[ 8] = -tmp;
91 
92  m[ 9] = 6.0 * (b3 - bias);
93  m[10] = 6.0 * bias;
94 
95  tmp = tension + 4.0 * (b2 + bias);
96  m[13] = tmp;
97  d = 1.0 / (tmp + 2.0 * b3 + 2.0);
98 
99  m[ 3] = m[14] = 2.0;
100  m[ 7] = m[11] = m[15] = 0.0;
101 
102  for (i=0; i < 16; i++) m[i] *= d;
103 }
104 
105 
106 static void
107 GetCardinal(fastf_t *m, const double tension)
108 {
109  m[ 1] = 2.0 - tension;
110  m[ 2] = tension - 2.0;
111  m[ 4] = 2.0 * tension;
112  m[ 5] = tension - 3.0;
113  m[ 6] = 3.0 - 2.0 * tension;
114  m[13] = 1.0;
115  m[ 3] = m[10] = tension;
116  m[ 0] = m[7] = m[ 8] = -tension;
117  m[ 9] = m[11] = m[12] = m[14] = m[15] = 0.0;
118 }
119 
120 
121 void
122 rt_dspline_matrix(fastf_t *m, const char *type, const double tension, const double bias)
123 
124 /* "Cardinal", "Catmull", "Beta" */
125 /* Cardinal tension of .5 is Catmull spline */
126 /* only for B spline */
127 {
128  if (!bu_strncasecmp(type, "Cardinal", 8)) {
129  GetCardinal(m, tension);
130  } else if (!bu_strncasecmp(type, "Catmull", 7)) {
131  GetCardinal(m, 0.5);
132  } else if (!bu_strncasecmp(type, "Beta", 4)) {
133  GetBeta(m, bias, tension);
134  } else {
135  bu_log("WARNING: %s:%d spline type \"%s\" unknown\n", __FILE__, __LINE__, type);
136  GetBeta(m, bias, tension);
137  }
138 }
139 
140 
141 double
142 rt_dspline4(fastf_t *m, double a, double b, double c, double d, double alpha)
143 /* spline matrix */
144 /* control pts */
145 /* point to interpolate at */
146 {
147  double p0, p1, p2, p3;
148 
149  p0 = m[ 0]*a + m[ 1]*b + m[ 2]*c + m[ 3]*d;
150  p1 = m[ 4]*a + m[ 5]*b + m[ 6]*c + m[ 7]*d;
151  p2 = m[ 8]*a + m[ 9]*b + m[10]*c + m[11]*d;
152  p3 = m[12]*a + m[13]*b + m[14]*c + m[15]*d;
153 
154  return p3 + alpha*(p2 + alpha*(p1 + alpha*p0));
155 }
156 
157 
158 void
159 rt_dspline4v(double *pt, const fastf_t *m, const double *a, const double *b, const double *c, const double *d, const int depth, const double alpha)
160 /* result */
161 /* spline matrix obtained with spline_matrix() */
162 /* knots */
163 
164 
165 /* number of values per knot */
166 /* 0 <= alpha <= 1 */
167 {
168  int i;
169  double p0, p1, p2, p3;
170 
171  for (i=0; i < depth; i++) {
172  p0 = m[ 0]*a[i] + m[ 1]*b[i] + m[ 2]*c[i] + m[ 3]*d[i];
173  p1 = m[ 4]*a[i] + m[ 5]*b[i] + m[ 6]*c[i] + m[ 7]*d[i];
174  p2 = m[ 8]*a[i] + m[ 9]*b[i] + m[10]*c[i] + m[11]*d[i];
175  p3 = m[12]*a[i] + m[13]*b[i] + m[14]*c[i] + m[15]*d[i];
176 
177  pt[i] = p3 + alpha*(p2 + alpha*(p1 + alpha*p0));
178  }
179 }
180 
181 void
182 rt_dspline_n(double *r, const fastf_t *m, const double *knots, const int nknots, const int depth, const double alpha)
183 /* result */
184 /* spline matrix */
185 /* knot values */
186 /* number of knots */
187 /* number of values per knot */
188 /* point on surface (0..1) to evaluate */
189 {
190  double *a, *b, *c, *d, x;
191  int nspans = nknots - 3;
192  int span;
193 
194  /* validate args */
195  if (nspans < 1 || depth < 1 || alpha < 0.0 || alpha > 1.0 ||
196  !r || !knots)
197  bu_bomb("invalid args given to rt_dspline_r");
198 
199 
200  /* compute which knots (span) we're going to interpolate */
201 
202 
203  x = alpha * nspans;
204  span = (int)x;
205  if (span >= nspans) span = nspans - 1;
206  x -= span;
207 
208  /* compute point (alpha 0..1) within this span */
209 
210  a = (double *)&knots[span*depth];
211  b = a+depth;
212  c = b+depth;
213  d = c+depth;
214 
215  rt_dspline4v(r, m, a, b, c, d, depth, x);
216 
217 }
218 
219 
220 /*
221  * Local Variables:
222  * mode: C
223  * tab-width: 8
224  * indent-tabs-mode: t
225  * c-file-style: "stroustrup"
226  * End:
227  * ex: shiftwidth=4 tabstop=8
228  */
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
Header file for the BRL-CAD common definitions.
int bu_strncasecmp(const char *string1, const char *string2, size_t n)
Definition: str.c:239
void rt_dspline_matrix(fastf_t *m, const char *type, const double tension, const double bias)
Definition: dspline.c:122
void rt_dspline4v(double *pt, const fastf_t *m, const double *a, const double *b, const double *c, const double *d, const int depth, const double alpha)
Definition: dspline.c:159
ustring alpha
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
double fastf_t
Definition: defines.h:300
void rt_dspline_n(double *r, const fastf_t *m, const double *knots, const int nknots, const int depth, const double alpha)
Definition: dspline.c:182
double rt_dspline4(fastf_t *m, double a, double b, double c, double d, double alpha)
Definition: dspline.c:142