BRL-CAD
oslo_calc.c
Go to the documentation of this file.
1 /* O S L O _ C A L C . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1986-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/oslo_calc.c
23  *
24  * Calculate the Oslo refinement matrix.
25  *
26  * This algorithm was taken from the paper:
27  *
28  * "Making the Oslo Algorithm More Efficient" by T. Lyche and K. Morken
29  *
30  * The algorithm referenced in the paper is algorithm 1 since we will
31  * be dealing mostly with surfaces. This routine computes the
32  * refinement matrix and returns a oslo structure which will allow a
33  * new curve or surface to be built.
34  *
35  * Since we only want the last row of the alpha's as outlined in the
36  * paper we can use a one dimensional array for the ah.
37  */
38 
39 #include "common.h"
40 
41 #include "bio.h"
42 
43 #include "vmath.h"
44 #include "raytrace.h"
45 #include "nurb.h"
46 
47 #define AMAX(i, j) ((i) > (j) ? (i) : (j))
48 #define AMIN(i, j) ((i) < (j) ? (i) : (j))
49 
50 struct oslo_mat *
51 rt_nurb_calc_oslo(register int order, register const struct knot_vector *tau_kv, register struct knot_vector *t_kv, struct resource *res)
52 
53 /* old knot vector */
54 /* new knot vector */
55 
56 {
57  register fastf_t *t_p;
58  register const fastf_t *tau_p;
59  fastf_t ah[20];
60  fastf_t newknots[20]; /* new knots */
61  register int j; /* d(j), j = 0 : # of new ctl points */
62  int mu, /* mu: tau[mu] <= t[j] < tau[mu+1]*/
63  muprim,
64  v, /* Nu value (order of matrix) */
65  p,
66  iu, /* upper bound loop counter */
67  il, /* lower bound loop counter */
68  ih,
69  n1; /* upper bound of t knot vector - order*/
70 
71  fastf_t tj;
72 
73  struct oslo_mat * head, * o_ptr, *new_o;
74 
75  if (res) RT_CK_RESOURCE(res);
76 
77  n1 = t_kv->k_size - order;
78 
79  t_p = t_kv->knots;
80  tau_p = tau_kv->knots;
81 
82  mu = 0; /* initialize mu */
83 
84  head = (struct oslo_mat *) bu_malloc (
85  sizeof(struct oslo_mat),
86  "rt_nurb_calc_oslo: oslo mat head");
87 
88  o_ptr = head;
89 
90  for (j = 0; j < n1; j++) {
91  register int i;
92 
93  if (j != 0) {
94  new_o = (struct oslo_mat *) bu_malloc (
95  sizeof(struct oslo_mat),
96  "rt_nurb_calc_oslo: oslo mat struct");
97 
98  o_ptr->next = new_o;
99  o_ptr = new_o;
100  }
101 
102  /* find the bounding mu */
103  while (tau_p[mu + 1] <= t_p[j]) mu++;
104 
105  muprim = mu;
106 
107  i = j + 1;
108 
109  while (ZERO(t_p[i] - tau_p[muprim]) && i < (j + order)) {
110  i++;
111  muprim--;
112  }
113 
114  ih = muprim + 1;
115 
116  for (v = 0, p = 1; p < order; p++) {
117  if (ZERO(t_p[j + p] - tau_p[ih]))
118  ih++;
119  else
120  newknots[++v - 1] = t_p[j + p];
121  }
122 
123  ah[order-1] = 1.0;
124 
125  for (p = 1; p <= v; p++) {
126 
127  fastf_t beta1;
128  int o_m;
129 
130  beta1 = 0.0;
131  o_m = order - muprim;
132 
133  tj = newknots[p-1];
134 
135  if (p > muprim) {
136  beta1 = ah[o_m];
137  beta1 = ((tj - tau_p[0]) * beta1) /
138  (tau_p[p + order - v] - tau_p[0]);
139  }
140  i = muprim - p + 1;
141  il = AMAX (1, i);
142  i = n1 - 1 + v - p;
143  iu = AMIN (muprim, i);
144 
145  for (i = il; i <= iu; i++) {
146  fastf_t d1, d2;
147  fastf_t beta;
148 
149  d1 = tj - tau_p[i];
150  d2 = tau_p[i + p + order - v - 1] - tj;
151 
152  beta = ah[i + o_m - 1] / (d1 + d2);
153 
154  ah[i + o_m - 2] = d2 * beta + beta1;
155  beta1 = d1 * beta;
156  }
157 
158  ah[iu + o_m - 1] = beta1;
159 
160  if (iu < muprim) {
161  register fastf_t kkk;
162  register fastf_t ahv;
163 
164  kkk = tau_p[n1 - 1 + order];
165  ahv = ah[iu + o_m];
166  ah[iu + o_m - 1] =
167  beta1 + (kkk - tj) *
168  ahv / (kkk - tau_p[iu + 1]);
169  }
170  }
171 
172  o_ptr->o_vec = (fastf_t *) bu_malloc (sizeof(fastf_t) * (v+1),
173  "rt_nurb_calc_oslo: oslo vector");
174 
175  o_ptr->offset = AMAX(muprim -v, 0);
176  o_ptr->osize = v;
177 
178  for (i = v, p = 0; i >= 0; i--)
179  o_ptr->o_vec[p++] = ah[(order-1) - i];
180  }
181 
182  o_ptr->next = (struct oslo_mat*) 0;
183  return head;
184 }
185 
186 
187 /**
188  * For debugging purposes only
189  */
190 void
191 rt_nurb_pr_oslo(struct oslo_mat *om)
192 {
193  struct oslo_mat * omp;
194  int j;
195 
196  for (omp = om; omp!= (struct oslo_mat *) 0; omp = omp->next) {
197  fprintf(stderr, "%p offset %d osize %d next %p\n",
198  (void *)omp, omp->offset, omp->osize,
199  (void *)omp->next);
200 
201  fprintf(stderr, "\t%f", omp->o_vec[0]);
202 
203  for (j = 1; j <= omp->osize; j++)
204  fprintf(stderr, "\t%f", omp->o_vec[j]);
205  fprintf(stderr, "\n");
206  }
207 }
208 
209 
210 /**
211  * Free up the structures and links for the oslo matrix.
212  */
213 void
214 rt_nurb_free_oslo(struct oslo_mat *om, struct resource *res)
215 {
216  register struct oslo_mat * omp;
217 
218  if (res) RT_CK_RESOURCE(res);
219 
220  while (om != (struct oslo_mat *) 0) {
221  omp = om;
222  om = om->next;
223  bu_free((char *)omp->o_vec, "rt_nurb_free_oslo: ovec");
224  bu_free((char *)omp, "rt_nurb_free_oslo: struct oslo");
225  }
226 }
227 
228 
229 /** @} */
230 /*
231  * Local Variables:
232  * mode: C
233  * tab-width: 8
234  * indent-tabs-mode: t
235  * c-file-style: "stroustrup"
236  * End:
237  * ex: shiftwidth=4 tabstop=8
238  */
Header file for the BRL-CAD common definitions.
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
void rt_nurb_free_oslo(struct oslo_mat *om, struct resource *res)
Definition: oslo_calc.c:214
#define ZERO(val)
Definition: units.c:38
#define RT_CK_RESOURCE(_p)
Definition: raytrace.h:1490
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
void rt_nurb_pr_oslo(struct oslo_mat *om)
Definition: oslo_calc.c:191
double fastf_t
Definition: defines.h:300
#define AMIN(i, j)
Definition: oslo_calc.c:48
struct oslo_mat * rt_nurb_calc_oslo(register int order, register const struct knot_vector *tau_kv, register struct knot_vector *t_kv, struct resource *res)
Definition: oslo_calc.c:51
#define AMAX(i, j)
Definition: oslo_calc.c:47