BRL-CAD
nurb_ray.c
Go to the documentation of this file.
1 /* N U R B _ R A Y . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1991-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_ray.c
23  *
24  * Functions which support the ray intersection for surfaces.
25  *
26  */
27 /** @} */
28 
29 #include "common.h"
30 
31 #include "bio.h"
32 
33 #include "vmath.h"
34 #include "nmg.h"
35 #include "raytrace.h"
36 #include "nurb.h"
37 
38 
39 void rt_nurb_pbound(struct face_g_snurb *srf, fastf_t *vmin, fastf_t *vmax);
40 
41 struct face_g_snurb *
42 rt_nurb_project_srf(const struct face_g_snurb *srf, fastf_t *plane1, fastf_t *plane2, struct resource *res)
43 {
44 
45  register struct face_g_snurb *psrf;
46  register fastf_t *mp1, *mp2;
47  int n_pt_type;
48  int rational;
49  int i;
50 
51  if (RTG.NMG_debug & DEBUG_RT_ISECT)
52  bu_log("rt_nurb_project_srf: projecting surface, planes = (%g %g %g %g) (%g %g %g %g)\n",
53  V4ARGS(plane1), V4ARGS(plane2));
54 
55  rational = RT_NURB_IS_PT_RATIONAL(srf->pt_type);
56 
57  n_pt_type = RT_NURB_MAKE_PT_TYPE(2, RT_NURB_PT_PROJ, 0);
58 
59  psrf = (struct face_g_snurb *) rt_nurb_new_snurb(srf->order[0], srf->order[1],
60  srf->u.k_size, srf->v.k_size,
61  srf->s_size[0], srf->s_size[1], n_pt_type, res);
62 
63  psrf->dir = RT_NURB_SPLIT_COL;
64 
65  for (i = 0; i < srf->u.k_size; i++) {
66  psrf->u.knots[i] = srf->u.knots[i];
67  }
68 
69  for (i = 0; i < srf->v.k_size; i++) {
70  psrf->v.knots[i] = srf->v.knots[i];
71  }
72 
73  mp1 = srf->ctl_points;
74  mp2 = psrf->ctl_points;
75 
76  for (i = 0; i < srf->s_size[0] * srf->s_size[1]; i++) {
77 
78  if (rational) {
79  mp2[0] = (mp1[0] / mp1[3] * plane1[0] +
80  mp1[1] / mp1[3] * plane1[1] +
81  mp1[2] / mp1[3] * plane1[2] - plane1[3]) *
82  mp1[3];
83  mp2[1] = (mp1[0] / mp1[3] * plane2[0] +
84  mp1[1] / mp1[3] * plane2[1] +
85  mp1[2] / mp1[3] * plane2[2] - plane2[3]) *
86  mp1[3];
87  } else {
88  mp2[0] = mp1[0] * plane1[0] + mp1[1] * plane1[1] +
89  mp1[2] * plane1[2] - plane1[3];
90  mp2[1] = mp1[0] * plane2[0] + mp1[1] * plane2[1] +
91  mp1[2] * plane2[2] - plane2[3];
92  }
93 
94  if (RTG.NMG_debug & DEBUG_RT_ISECT) {
95  if (rational)
96  bu_log("\tmesh pt (%g %g %g %g), becomes (%g %g)\n", V4ARGS(mp1), mp2[0], mp2[1]);
97  else
98  bu_log("\tmesh pt (%g %g %g), becomes (%g %g)\n", V3ARGS(mp1), mp2[0], mp2[1]);
99  }
100 
101  mp1 += RT_NURB_EXTRACT_COORDS(srf->pt_type);
102  mp2 += RT_NURB_EXTRACT_COORDS(psrf->pt_type);
103  }
104 
105  return (struct face_g_snurb *) psrf;
106 }
107 
108 
109 /**
110  * This routine should go away and be changed into a macro
111  * but for now I want to be able to useit with debugging.
112  * - Paul
113  */
114 #define FINDZERO(x0, x1, y0, y1) (x0 - y0 * (x1 - x0) / (y1-y0))
115 
118 };
119 
120 
124 };
125 
126 #if !defined(SIGN)
127 #define SIGN(a) ((a < 0.0)? -1 : 1)
128 #endif
129 
130 void
131 rt_nurb_clip_srf(const struct face_g_snurb *srf, int dir, fastf_t *min, fastf_t *max)
132 {
133  struct internal_convex_hull ch[20]; /* max order is 10 */
134  register fastf_t * mp1;
135  fastf_t * p1, *p2, *p3, *p4; /* corner points of the mesh */
136  fastf_t v1[2], v2[2], v3[2]; /* vectors from corners */
137  struct internal_line l1;
138  fastf_t norm;
139  fastf_t value;
140  int i;
141  register int j;
142  int k;
143  int coords;
144  int col_size, row_size;
145 
146  col_size = srf->s_size[1];
147  row_size = srf->s_size[0];
148 
149  coords = RT_NURB_EXTRACT_COORDS(srf->pt_type);
150 
151  p1 = srf->ctl_points;
152  p2 = srf->ctl_points + coords * (col_size - 1);
153  p3 = srf->ctl_points + (coords * col_size *
154  (row_size - 1));
155  p4 = srf->ctl_points + (coords * col_size *
156  (row_size - 1)) +
157  ((col_size - 1) * coords);
158 
159  if (dir == RT_NURB_SPLIT_ROW) {
160  v1[0] = p1[0] - p3[0];
161  v1[1] = p1[1] - p3[1];
162 
163  v2[0] = p2[0] - p4[0];
164  v2[1] = p2[1] - p4[1];
165  } else {
166  v1[0] = p1[0] - p2[0];
167  v1[1] = p1[1] - p2[1];
168 
169  v2[0] = p3[0] - p4[0];
170  v2[1] = p3[1] - p4[1];
171  }
172 
173  v3[0] = v1[0] + v2[0];
174  v3[1] = v1[1] + v1[1];
175 
176  norm = sqrt(v3[1] * v3[1] + v3[0] * v3[0]);
177  l1.a = v3[1] / norm;
178  l1.b = -v3[0] / norm;
179 
180  *min = 1.0e8;
181  *max = -1.0e8;
182 
183  if (dir == RT_NURB_SPLIT_ROW) {
184  for (i = 0; i < col_size; i++) {
185  ch[i].param = (fastf_t) i / (col_size - 1.0);
186  ch[i].min = 1.0e8;
187  ch[i].max = -1.0e8;
188  }
189 
190  mp1 = srf->ctl_points;
191 
192  for (i = 0; i < row_size; i++) {
193  for (j = 0; j < col_size; j++) {
194  value = - (mp1[0] * l1.a + mp1[1] * l1.b);
195  if (value <= ch[j].min)
196  ch[j].min = value;
197  if (value >= ch[j].max)
198  ch[j].max = value;
199  mp1 += coords;
200  }
201  }
202 
203  for (k = 0; k < col_size - 1; k++)
204  for (j = k+1; j < col_size; j++) {
205  fastf_t d;
206  fastf_t param1, param2;
207 
208  param1 = ch[k].param;
209  param2 = ch[j].param;
210 
211  d = FINDZERO(param1, param2, ch[k].max, ch[j].max);
212  if (d <= *min) *min = d * .99;
213  if (d >= *max) *max = d * .99 + .01;
214 
215  d = FINDZERO(param1, param2, ch[k].min, ch[j].min);
216  if (d <= *min) *min = d * .99;
217  if (d >= *max) *max = d * .99 + .01;
218  }
219 
220  if (*min <= 0.0)
221  *min = 0.0;
222  if (*max >= 1.0)
223  *max = 1.0;
224  if (SIGN(ch[0].min) != SIGN(ch[0].max))
225  *min = 0.0;
226  i = SIGN(ch[col_size -1].min);
227  j = SIGN(ch[col_size -1].max);
228  if (i != j)
229  *max = 1.0;
230  } else {
231  for (i = 0; i < row_size; i++) {
232  ch[i].param = (fastf_t) i / (row_size - 1.0);
233  ch[i].min = 1.0e8;
234  ch[i].max = -1.0e8;
235  }
236 
237 
238  for (i = 0; i < col_size; i++) {
239  int stride;
240 
241  stride = coords * col_size;
242 
243  mp1 = srf->ctl_points + i * coords;
244  for (j = 0; j < row_size; j++) {
245  value = - (mp1[0] * l1.a + mp1[1] * l1.b);
246  if (value <= ch[j].min)
247  ch[j].min = value;
248  if (value >= ch[j].max)
249  ch[j].max = value;
250  mp1 += stride;
251  }
252  }
253 
254  for (k = 0; k < row_size - 1; k++)
255  for (j = k+1; j < row_size; j++) {
256  fastf_t d;
257  fastf_t param1, param2;
258 
259  param1 = ch[k].param;
260  param2 = ch[j].param;
261 
262  d = FINDZERO(param1, param2, ch[k].max, ch[j].max);
263  if (d <= *min) *min = d * .99;
264  if (d >= *max) *max = d * .99 + .01;
265 
266  d = FINDZERO(param1, param2, ch[k].min, ch[j].min);
267  if (d <= *min) *min = d * .99;
268  if (d >= *max) *max = d * .99 + .01;
269  }
270  if (*min <= 0.0)
271  *min = 0.0;
272  if (*max >= 1.0)
273  *max = 1.0;
274  if (SIGN(ch[0].min) != SIGN(ch[0].max))
275  *min = 0.0;
276  i = SIGN(ch[row_size-1 ].min);
277  j = SIGN(ch[row_size -1].max);
278  if (i != j)
279  *max = 1.0; }
280 }
281 
282 
283 struct face_g_snurb *
284 rt_nurb_region_from_srf(const struct face_g_snurb *srf, int dir, fastf_t param1, fastf_t param2, struct resource *res)
285 {
286  register int i;
287  struct face_g_snurb *region;
288  struct knot_vector new_knots;
289 
290  fastf_t *knot_vec = NULL;
291  size_t maxorder = FMAX(srf->order[0], srf->order[1]);
292  knot_vec = (fastf_t *)bu_calloc(maxorder * 2, sizeof(fastf_t), "knot vector");
293 
294  /* Build the new knot vector in a local array, which gets copied
295  * later in rt_nurb_s_refine(). */
296  new_knots.knots = &knot_vec[0];
297 
298  if (dir == RT_NURB_SPLIT_ROW) {
299  new_knots.k_size = srf->order[0] * 2;
300 
301  for (i = 0; i < srf->order[0]; i++) {
302  knot_vec[i] = param1;
303  knot_vec[i+srf->order[0]] = param2;
304  }
305  } else {
306  new_knots.k_size = srf->order[1] * 2;
307 
308  for (i = 0; i < srf->order[1]; i++) {
309  knot_vec[i] = param1;
310  knot_vec[i+srf->order[1]] = param2;
311  }
312  }
313 
314  region = rt_nurb_s_refine(srf, dir, &new_knots, res);
315  bu_free(knot_vec, "knot vector");
316 
317  return region;
318 }
319 
320 
321 struct rt_nurb_uv_hit *
322 rt_nurb_intersect(const struct face_g_snurb *srf, fastf_t *plane1, fastf_t *plane2, double uv_tol, struct resource *res, struct bu_list *plist)
323 {
324  struct rt_nurb_uv_hit * h;
325  struct face_g_snurb * psrf,
326  * osrf;
327  int dir,
328  sub;
329 
330  point_t vmin,
331  vmax;
332  fastf_t u[2],
333  v[2];
334  struct bu_list rni_plist;
335 
336  NMG_CK_SNURB(srf);
337 
338  h = (struct rt_nurb_uv_hit *) 0;
339  if (plist == NULL) {
340  plist = &rni_plist;
341  BU_LIST_INIT(plist);
342  }
343 
344  /* project the surface to a 2 dimensional problem */
345  /* NOTE that this gives a single snurb back, NOT a list */
346  psrf = rt_nurb_project_srf(srf, plane2, plane1, res);
347  psrf->dir = 1;
348  BU_LIST_APPEND(plist, &psrf->l);
349 
350  if (RT_G_DEBUG & DEBUG_SPLINE)
351  rt_nurb_s_print("srf", psrf);
352 
353  /* This list starts out with only a single snurb, but more may be
354  * added on as work progresses.
355  */
356  while (BU_LIST_WHILE(psrf, face_g_snurb, plist)) {
357  int flat;
358 
359  BU_LIST_DEQUEUE(&psrf->l);
360  NMG_CK_SNURB(psrf);
361  sub = 0;
362  flat = 0;
363  dir = psrf->dir;
364 
365  while (!flat) {
366  fastf_t smin = 0.0, smax = 0.0;
367 
368  sub++;
369  dir = (dir == 0)?1:0; /* change direction */
370 
371  if (RT_G_DEBUG & DEBUG_SPLINE)
372  rt_nurb_s_print("psrf", psrf);
373 
374  rt_nurb_pbound(psrf, vmin, vmax);
375 
376  /* Check for origin to be included in the bounding box */
377  if (!(vmin[0] <= 0.0 && vmin[1] <= 0.0 &&
378  vmax[0] >= 0.0 && vmax[1] >= 0.0)) {
379  if (RT_G_DEBUG & DEBUG_SPLINE)
380  bu_log("this srf doesn't include the origin\n");
381  flat = 1;
382  rt_nurb_free_snurb(psrf, res);
383  continue;
384  }
385 
386  rt_nurb_clip_srf(psrf, dir, &smin, &smax);
387 
388  if ((smax - smin) > .8) {
389  struct rt_nurb_uv_hit *hp;
390 
391  /* Split surf, requeue both sub-surfs at head */
392  /* New surfs will have same dir as arg, here */
393  if (RT_G_DEBUG & DEBUG_SPLINE)
394  bu_log("splitting this surface\n");
395  rt_nurb_s_split(plist, psrf, dir, res);
396  rt_nurb_free_snurb(psrf, res);
397 
398  hp = rt_nurb_intersect(srf, plane1, plane2, uv_tol, res, plist);
399  return hp;
400  }
401  if (smin > 1.0 || smax < 0.0) {
402  if (RT_G_DEBUG & DEBUG_SPLINE)
403  bu_log("eliminating this surface (smin=%g, smax=%g)\n", smin, smax);
404  flat = 1;
405  rt_nurb_free_snurb(psrf, res);
406  continue;
407  }
408  if (dir == RT_NURB_SPLIT_ROW) {
409  smin = (1.0 - smin) * psrf->u.knots[0] +
410  smin * psrf->u.knots[
411  psrf->u.k_size -1];
412  smax = (1.0 - smax) * psrf->u.knots[0] +
413  smax * psrf->u.knots[
414  psrf->u.k_size -1];
415  } else {
416  smin = (1.0 - smin) * psrf->v.knots[0] +
417  smin * psrf->v.knots[
418  psrf->v.k_size -1];
419  smax = (1.0 - smax) * psrf->v.knots[0] +
420  smax * psrf->v.knots[
421  psrf->v.k_size -1];
422  }
423 
424  osrf = psrf;
425  psrf = (struct face_g_snurb *) rt_nurb_region_from_srf(
426  osrf, dir, smin, smax, res);
427 
428  psrf->dir = dir;
429  rt_nurb_free_snurb(osrf, res);
430 
431  if (RT_G_DEBUG & DEBUG_SPLINE) {
432  bu_log("After call to rt_nurb_region_from_srf() (smin=%g, smax=%g)\n", smin, smax);
433  rt_nurb_s_print("psrf", psrf);
434  }
435 
436  u[0] = psrf->u.knots[0];
437  u[1] = psrf->u.knots[psrf->u.k_size -1];
438 
439  v[0] = psrf->v.knots[0];
440  v[1] = psrf->v.knots[psrf->v.k_size -1];
441 
442  if ((u[1] - u[0]) < uv_tol && (v[1] - v[0]) < uv_tol) {
443  struct rt_nurb_uv_hit * hit;
444 
445  if (RT_G_DEBUG & DEBUG_SPLINE) {
446  fastf_t p1[4], p2[4];
447  int coords;
448  vect_t diff;
449 
450  coords = RT_NURB_EXTRACT_COORDS(srf->pt_type);
451  rt_nurb_s_eval(srf, u[0], v[0], p1);
452  rt_nurb_s_eval(srf, u[1], v[1], p2);
453 
454  if (RT_NURB_IS_PT_RATIONAL(srf->pt_type)) {
455  fastf_t inv_w;
456 
457  inv_w = 1.0 / p1[coords-1];
458  VSCALE(p1, p1, inv_w);
459 
460  inv_w = 1.0 / p2[coords-1];
461  VSCALE(p2, p2, inv_w);
462  }
463 
464  VSUB2(diff, p1, p2);
465  bu_log("Precision of hit point = %g (%f %f %f) <-> (%f %f %f)\n",
466  MAGNITUDE(diff), V3ARGS(p1), V3ARGS(p2));
467  }
468 
469  hit = (struct rt_nurb_uv_hit *) bu_malloc(
470  sizeof(struct rt_nurb_uv_hit), "hit");
471 
472  hit->next = (struct rt_nurb_uv_hit *)0;
473  hit->sub = sub;
474  hit->u = (u[0] + u[1])/2.0;
475  hit->v = (v[0] + v[1])/2.0;
476 
477  if (h == (struct rt_nurb_uv_hit *)0)
478  h = hit;
479  else {
480  hit->next = h;
481  h = hit;
482  }
483  flat = 1;
484  rt_nurb_free_snurb(psrf, res);
485  }
486  if ((u[1] - u[0]) > (v[1] - v[0]))
487  dir = 1;
488  else dir = 0;
489  }
490  }
491 
492  return (struct rt_nurb_uv_hit *)h;
493 }
494 
495 
496 void
497 rt_nurb_pbound(struct face_g_snurb *srf, fastf_t *vmin, fastf_t *vmax)
498 {
499  register fastf_t * ptr;
500  register int coords;
501  int i;
502 
503  vmin[0] = vmin[1] = vmin[2] = INFINITY;
504  vmax[0] = vmax[1] = vmax[2] = -INFINITY;
505 
506  ptr = srf->ctl_points;
507 
508  coords = RT_NURB_EXTRACT_COORDS(srf->pt_type);
509 
510  for (i = (srf->s_size[RT_NURB_SPLIT_ROW] *
511  srf->s_size[RT_NURB_SPLIT_COL]); i > 0; i--) {
512  V_MIN((vmin[0]), (ptr[0]));
513  V_MAX((vmax[0]), (ptr[0]));
514 
515  V_MIN((vmin[1]), (ptr[1]));
516  V_MAX((vmax[1]), (ptr[1]));
517 
518  ptr += coords;
519  }
520 }
521 
522 
523 /*
524  * Local Variables:
525  * mode: C
526  * tab-width: 8
527  * indent-tabs-mode: t
528  * c-file-style: "stroustrup"
529  * End:
530  * ex: shiftwidth=4 tabstop=8
531  */
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
Definition: list.h:118
struct face_g_snurb * rt_nurb_s_refine(const struct face_g_snurb *srf, int dir, struct knot_vector *kv, struct resource *res)
Definition: nurb_refine.c:45
struct face_g_snurb * rt_nurb_new_snurb(int u_order, int v_order, int n_u, int n_v, int n_rows, int n_cols, int pt_type, struct resource *res)
Definition: nurb_util.c:42
Header file for the BRL-CAD common definitions.
#define BU_LIST_APPEND(old, new)
Definition: list.h:197
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
#define RT_G_DEBUG
Definition: raytrace.h:1718
uint32_t NMG_debug
debug bits for NMG's see nmg.h
Definition: raytrace.h:1699
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
#define FINDZERO(x0, x1, y0, y1)
Definition: nurb_ray.c:114
#define V3ARGS(a)
Definition: color.c:56
fastf_t b
Definition: nurb_ray.c:117
#define SIGN(a)
Definition: nurb_ray.c:127
void rt_nurb_pbound(struct face_g_snurb *srf, fastf_t *vmin, fastf_t *vmax)
Definition: nurb_ray.c:497
void rt_nurb_s_print(char *c, const struct face_g_snurb *srf)
Definition: nurb_util.c:214
#define BU_LIST_WHILE(p, structure, hp)
Definition: list.h:410
fastf_t a
Definition: nurb_ray.c:117
#define BU_LIST_INIT(_hp)
Definition: list.h:148
struct face_g_snurb * rt_nurb_project_srf(const struct face_g_snurb *srf, fastf_t *plane1, fastf_t *plane2, struct resource *res)
Definition: nurb_ray.c:42
int hit(struct application *ap, struct partition *PartHeadp, struct seg *segs)
Definition: gqa.c:963
void rt_nurb_clip_srf(const struct face_g_snurb *srf, int dir, fastf_t *min, fastf_t *max)
Definition: nurb_ray.c:131
#define FMAX(a, b)
Definition: common.h:102
struct face_g_snurb * rt_nurb_region_from_srf(const struct face_g_snurb *srf, int dir, fastf_t param1, fastf_t param2, struct resource *res)
Definition: nurb_ray.c:284
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
void rt_nurb_s_split(struct bu_list *split_hd, const struct face_g_snurb *srf, int dir, struct resource *res)
Definition: nurb_split.c:56
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
#define DEBUG_SPLINE
9 Splines
Definition: raytrace.h:92
#define BU_LIST_DEQUEUE(cur)
Definition: list.h:209
struct rt_nurb_uv_hit * rt_nurb_intersect(const struct face_g_snurb *srf, fastf_t *plane1, fastf_t *plane2, double uv_tol, struct resource *res, struct bu_list *plist)
Definition: nurb_ray.c:322
double fastf_t
Definition: defines.h:300
struct rt_g RTG
Definition: globals.c:39