BRL-CAD
nurb_interp.c
Go to the documentation of this file.
1 /* N U R B _ I N T E R P . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1994-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_interp.c
23  *
24  * Interpolation routines for fitting NURB curves and and surfaces to
25  * existing data.
26  *
27  */
28 /** @} */
29 
30 #include "common.h"
31 
32 #include <string.h>
33 #include "bio.h"
34 
35 #include "vmath.h"
36 #include "raytrace.h"
37 #include "nurb.h"
38 
39 
40 void
41 rt_nurb_nodes(fastf_t *nodes, const struct knot_vector *knots, int order)
42 {
43  int i, j;
44  fastf_t sum;
45 
46  for (i = 0; i < knots->k_size -order; i++) {
47 
48  sum = 0.0;
49 
50  for (j = 1; j <= order -1; j++) {
51  sum += knots->knots[i+j];
52  }
53  nodes[i] = sum/(order -1);
54  }
55 }
56 
57 
58 void
59 rt_nurb_interp_mat(fastf_t *imat, struct knot_vector *knots, fastf_t *nodes, int order, int dim)
60 {
61  int i, j;
62  int ptr;
63 
64  ptr = 0;
65 
66  for (i = 0; i < dim; i++)
67  for (j = 0; j < dim; j++) {
68  imat[ptr] = rt_nurb_basis_eval(knots, j, order, nodes[i]);
69  ptr++;
70  }
71 
72  imat[ptr-1] = 1.0;
73 }
74 
75 
76 /**
77  * main routine for interpolation of curves
78  */
79 void
80 rt_nurb_cinterp(struct edge_g_cnurb *crv, int order, const fastf_t *data, int n)
81 {
82  fastf_t * interp_mat;
83  fastf_t * nodes;
84  fastf_t *local_data;
85 
86  /* Create Data memory and fill in curve structs */
87 
88  interp_mat = (fastf_t *) bu_malloc(n * n * sizeof(fastf_t),
89  "rt_nurb_interp: interp_mat");
90 
91  nodes = (fastf_t *) bu_malloc(n * sizeof(fastf_t), "rt_nurb_interp:nodes");
92  local_data = (fastf_t *)bu_malloc(n * 3 * sizeof(fastf_t), "rt_nurb_interp() local_data[]");
93 
94  crv->ctl_points = (fastf_t *) bu_malloc(n * 3 * sizeof(fastf_t),
95  "solution");
96 
97  crv->order = order;
98  crv->c_size = n;
99  crv->pt_type = RT_NURB_MAKE_PT_TYPE(3, RT_NURB_PT_XYZ, 0);
100 
101  /* First set up Curve data structs */
102  /* For now we will assume that all parameterizations are uniform */
103 
104  rt_nurb_kvknot(&crv->k, order, 0.0, 1.0, (n - order), (struct resource *)NULL);
105 
106  /* Calculate Nodes at which the data points will be evaluated in
107  * the curve
108  */
109 
110  rt_nurb_nodes(nodes, &crv->k, order);
111 
112  /* use the node values to create the interpolation matrix which is
113  * a diagonal matrix
114  */
115 
116  rt_nurb_interp_mat(interp_mat, &crv->k, nodes, order, n);
117 
118  /* Solve the system of equations to get the control points Because
119  * rt_nurb_solve needs to modify the data as it works, and it
120  * wouldn't be polite to trash our caller's data, make a local
121  * copy. This creates the final ctl_points[] array.
122  */
123  memcpy((char *)local_data, (char *)data, n * 3 * sizeof(fastf_t));
124  rt_nurb_solve(interp_mat, local_data, crv->ctl_points, n, 3);
125 
126  /* Free up node and interp_mat storage */
127 
128  bu_free((char *) interp_mat, "rt_nurb_cinterp: interp_mat");
129  bu_free((char *) nodes, "rt_nurb_cinterp: nodes");
130  bu_free((char *) local_data, "rt_nurb_cinterp() local_data[]");
131 
132  /* All done, The resulting crv now interpolates the data */
133 }
134 
135 
136 /**
137  * Interpolate the 2-D grid of data values and fit a B-spline surface
138  * to it.
139  *
140  * This is done in two steps:
141  *
142  * 1) Fit a curve to the data in each row.
143  *
144  * 2) Fit a curve to the control points from step 1 in each column.
145  * The result is a mesh of control points which defines the surface.
146  *
147  * Input data is assumed to be a 3-tuple of (X, Y, Z) where Z is the
148  * independent variable being interpolated to make the surface.
149  */
150 void
151 rt_nurb_sinterp(struct face_g_snurb *srf, int order, const fastf_t *data, int ymax, int xmax)
152 
153 
154  /* data[x, y] */
155  /* nrow = max Y */
156  /* ncol = max X */
157 {
158  int x;
159  int y;
160  struct edge_g_cnurb *crv; /* array of cnurbs */
161  fastf_t *tmp;
162  fastf_t *cpt; /* surface control point pointer */
163 
164  /* Build the resultant surface structure */
165  srf->order[0] = srf->order[1] = order;
166  srf->dir = 0;
167  srf->s_size[0] = xmax;
168  srf->s_size[1] = ymax;
169  srf->l.magic = NMG_FACE_G_SNURB_MAGIC;
170  srf->pt_type = RT_NURB_MAKE_PT_TYPE(3, RT_NURB_PT_XYZ, RT_NURB_PT_NONRAT);
171 
172  /* the U knot vector replates to the points in a row therefore you
173  * want to determine how many cols there are similar for the V knot
174  * vector
175  */
176 
177  rt_nurb_kvknot(&srf->u, order, 0.0, 1.0, ymax - order, (struct resource *)NULL);
178  rt_nurb_kvknot(&srf->v, order, 0.0, 1.0, xmax - order, (struct resource *)NULL);
179 
180  srf->ctl_points = (fastf_t *) bu_malloc(
181  sizeof(fastf_t) * xmax * ymax * 3,
182  "rt_nurb_sinterp() surface ctl_points[]");
183  cpt = &srf->ctl_points[0];
184 
185 /* _col is X, _row is Y */
186 #define NVAL(_col, _row) data[((_row)*xmax+(_col))*3]
187 
188  crv = (struct edge_g_cnurb *)bu_calloc(sizeof(struct edge_g_cnurb), ymax,
189  "rt_nurb_sinterp() crv[]");
190 
191  /* Interpolate the data across the rows, fitting a curve to each. */
192  for (y = 0; y < ymax; y++) {
193  crv[y].l.magic = NMG_EDGE_G_CNURB_MAGIC;
194  /* Build curve from from (0, y) to (xmax-1, y) */
195  rt_nurb_cinterp(&crv[y], order, &NVAL(0, y), xmax);
196  }
197 #undef NVAL
198 
199  tmp = (fastf_t *)bu_malloc(sizeof(fastf_t)*3 * ymax,
200  "rt_nurb_sinterp() tmp[]");
201  for (x = 0; x < xmax; x++) {
202  struct edge_g_cnurb ncrv;
203 
204  /* copy the curve ctl points into col major format */
205  for (y = 0; y < ymax; y++) {
206  VMOVE(&tmp[y*3], &crv[y].ctl_points[x*3]);
207  }
208 
209  /* Interpolate the curve interpolates, giving rows of a surface */
210  ncrv.l.magic = NMG_EDGE_G_CNURB_MAGIC;
211  rt_nurb_cinterp(&ncrv, order, tmp, ymax);
212 
213  /* Move new curve interpolations into snurb ctl_points[] */
214  for (y = 0; y < ymax*3; y++) {
215  *cpt++ = ncrv.ctl_points[y];
216  }
217  rt_nurb_clean_cnurb(&ncrv);
218  }
219  for (y = 0; y < ymax; y++) {
220  rt_nurb_clean_cnurb(&crv[y]);
221  }
222  bu_free((char *)crv, "crv[]");
223  bu_free((char *)tmp, "tmp[]");
224 }
225 
226 
227 /*
228  * Local Variables:
229  * mode: C
230  * tab-width: 8
231  * indent-tabs-mode: t
232  * c-file-style: "stroustrup"
233  * End:
234  * ex: shiftwidth=4 tabstop=8
235  */
void rt_nurb_nodes(fastf_t *nodes, const struct knot_vector *knots, int order)
Definition: nurb_interp.c:41
void rt_nurb_cinterp(struct edge_g_cnurb *crv, int order, const fastf_t *data, int n)
Definition: nurb_interp.c:80
#define NMG_FACE_G_SNURB_MAGIC
Definition: magic.h:126
Header file for the BRL-CAD common definitions.
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
COMPLEX data[64]
Definition: fftest.c:34
void rt_nurb_sinterp(struct face_g_snurb *srf, int order, const fastf_t *data, int ymax, int xmax)
Definition: nurb_interp.c:151
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
#define NMG_EDGE_G_CNURB_MAGIC
Definition: magic.h:121
void rt_nurb_interp_mat(fastf_t *imat, struct knot_vector *knots, fastf_t *nodes, int order, int dim)
Definition: nurb_interp.c:59
#define NVAL(_col, _row)
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
void rt_nurb_clean_cnurb(struct edge_g_cnurb *crv)
Definition: nurb_util.c:152
double fastf_t
Definition: defines.h:300
void rt_nurb_solve(fastf_t *mat_1, fastf_t *mat_2, fastf_t *solution, int dim, int coords)
Definition: nurb_solve.c:58
void rt_nurb_kvknot(register struct knot_vector *new_knots, int order, fastf_t lower, fastf_t upper, int num, struct resource *res)
Definition: nurb_knot.c:47
fastf_t rt_nurb_basis_eval(register struct knot_vector *knts, int interval, int order, fastf_t mu)
Definition: nurb_basis.c:63