BRL-CAD
libbrep_brep_tools.cpp
Go to the documentation of this file.
1 /* L I B B R E P _ B R E P _ T O O L S . C P P
2  * BRL-CAD
3  *
4  * Copyright (c) 2013-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 /** @file libbrep_brep_tools.cpp
21  *
22  * Utility routines for working with geometry.
23  */
24 
25 #include "common.h"
26 
27 #include <vector>
28 #include <iostream>
29 
30 #include "opennurbs.h"
31 #include "bu/log.h"
32 #include "libbrep_brep_tools.h"
33 
34 bool ON_NearZero(double val, double epsilon) {
35  return (val > -epsilon) && (val < epsilon);
36 }
37 
38 double ON_Curve_Get_Tangent(int direction, const ON_Curve* curve, double min, double max, double zero_tol) {
39  double mid;
40  bool tanmin;
41  ON_3dVector tangent;
42 
43  // first check end points
44  tangent = curve->TangentAt(max);
45  if (direction == 1 && ON_NearZero(tangent.y, zero_tol)) return max;
46  if (direction == 0 && ON_NearZero(tangent.x, zero_tol)) return max;
47  tangent = curve->TangentAt(min);
48  if (direction == 1 && ON_NearZero(tangent.y, zero_tol)) return min;
49  if (direction == 0 && ON_NearZero(tangent.y, zero_tol)) return min;
50 
51  tanmin = (tangent[direction] < 0.0);
52  while ((max-min) > zero_tol) {
53  mid = (max + min)/2.0;
54  tangent = curve->TangentAt(mid);
55  if (ON_NearZero(tangent[direction], zero_tol)) {
56  return mid;
57  }
58  if ((tangent[direction] < 0.0) == tanmin) {
59  min = mid;
60  } else {
61  max = mid;
62  }
63  }
64  return min;
65 }
66 
67 double ON_Curve_Get_Horizontal_Tangent(const ON_Curve* curve, double min, double max, double zero_tol) {
68  return ON_Curve_Get_Tangent(1, curve, min, max, zero_tol);
69 }
70 
71 double ON_Curve_Get_Vertical_Tangent(const ON_Curve* curve, double min, double max, double zero_tol) {
72  return ON_Curve_Get_Tangent(0, curve, min, max, zero_tol);
73 }
74 
75 int ON_Curve_Has_Tangent(const ON_Curve* curve, double ct_min, double ct_max, double t_tol) {
76 
77  bool tanx1, tanx2, x_changed;
78  bool tany1, tany2, y_changed;
79  bool slopex, slopey;
80  double xdelta, ydelta;
81  ON_3dVector tangent1, tangent2;
82  ON_3dPoint p1, p2;
83  ON_Interval t(ct_min, ct_max);
84 
85  tangent1 = curve->TangentAt(t[0]);
86  tangent2 = curve->TangentAt(t[1]);
87 
88  tanx1 = (tangent1[0] < 0.0);
89  tanx2 = (tangent2[0] < 0.0);
90  tany1 = (tangent1[1] < 0.0);
91  tany2 = (tangent2[1] < 0.0);
92 
93  x_changed =(tanx1 != tanx2);
94  y_changed =(tany1 != tany2);
95 
96  if (x_changed && y_changed) return 3; //horz & vert
97  if (x_changed) return 1;//need to get vertical tangent
98  if (y_changed) return 2;//need to find horizontal tangent
99 
100  p1 = curve->PointAt(t[0]);
101  p2 = curve->PointAt(t[1]);
102 
103  xdelta = (p2[0] - p1[0]);
104  slopex = (xdelta < 0.0);
105  ydelta = (p2[1] - p1[1]);
106  slopey = (ydelta < 0.0);
107 
108  // If we have no slope change
109  // in x or y, we have a tangent line
110  if (ON_NearZero(xdelta, t_tol) || ON_NearZero(ydelta, t_tol)) return 0;
111 
112  if ((slopex != tanx1) || (slopey != tany1)) return 3;
113 
114  return 0;
115 }
116 
117 bool ON_Surface_IsFlat(ON_Plane *frames, double f_tol)
118 {
119  double Ndot=1.0;
120 
121  for (int i=0; i<8; i++) {
122  for ( int j=i+1; j<9; j++) {
123  if ((Ndot = Ndot * frames[i].zaxis * frames[j].zaxis) < f_tol) {
124  return false;
125  }
126  }
127  }
128 
129  return true;
130 }
131 
132 bool ON_Surface_IsFlat_U(ON_Plane *frames, double f_tol)
133 {
134  // check surface normals in U direction
135  double Ndot = 1.0;
136  if ((Ndot=frames[0].zaxis * frames[1].zaxis) < f_tol) {
137  return false;
138  } else if ((Ndot=Ndot * frames[2].zaxis * frames[3].zaxis) < f_tol) {
139  return false;
140  } else if ((Ndot=Ndot * frames[5].zaxis * frames[7].zaxis) < f_tol) {
141  return false;
142  } else if ((Ndot=Ndot * frames[6].zaxis * frames[8].zaxis) < f_tol) {
143  return false;
144  }
145 
146  // check for U twist within plane
147  double Xdot = 1.0;
148  if ((Xdot=frames[0].xaxis * frames[1].xaxis) < f_tol) {
149  return false;
150  } else if ((Xdot=Xdot * frames[2].xaxis * frames[3].xaxis) < f_tol) {
151  return false;
152  } else if ((Xdot=Xdot * frames[5].xaxis * frames[7].xaxis) < f_tol) {
153  return false;
154  } else if ((Xdot=Xdot * frames[6].xaxis * frames[8].xaxis) < f_tol) {
155  return false;
156  }
157 
158  return true;
159 }
160 
161 bool ON_Surface_IsFlat_V(ON_Plane *frames, double f_tol)
162 {
163  // check surface normals in V direction
164  double Ndot = 1.0;
165  if ((Ndot=frames[0].zaxis * frames[3].zaxis) < f_tol) {
166  return false;
167  } else if ((Ndot=Ndot * frames[1].zaxis * frames[2].zaxis) < f_tol) {
168  return false;
169  } else if ((Ndot=Ndot * frames[5].zaxis * frames[6].zaxis) < f_tol) {
170  return false;
171  } else if ((Ndot=Ndot * frames[7].zaxis * frames[8].zaxis) < f_tol) {
172  return false;
173  }
174 
175  // check for V twist within plane
176  double Xdot = 1.0;
177  if ((Xdot=frames[0].xaxis * frames[3].xaxis) < f_tol) {
178  return false;
179  } else if ((Xdot=Xdot * frames[1].xaxis * frames[2].xaxis) < f_tol) {
180  return false;
181  } else if ((Xdot=Xdot * frames[5].xaxis * frames[6].xaxis) < f_tol) {
182  return false;
183  } else if ((Xdot=Xdot * frames[7].xaxis * frames[8].xaxis) < f_tol) {
184  return false;
185  }
186 
187  return true;
188 }
189 
190 bool ON_Surface_IsStraight(ON_Plane *frames, double s_tol)
191 {
192  double Xdot=1.0;
193 
194  for (int i=0; i<8; i++) {
195  for ( int j=i+1; j<9; j++) {
196  if ((Xdot = Xdot * frames[0].xaxis * frames[1].xaxis) < s_tol) {
197  return false;
198  }
199  }
200  }
201 
202  return true;
203 }
204 
205 /**
206  \brief Create surfaces and store their pointers in the t* arguments.
207 
208  For any pre-existing surface passed as one of the t* args, this is a no-op.
209 
210  @param t1 Pointer to pointer addressing first surface
211  @param t2 Pointer to pointer addressing second surface
212  @param t3 Pointer to pointer addressing third surface
213  @param t4 Pointer to pointer addressing fourth surface
214 */
216  ON_Surface **t1,
217  ON_Surface **t2,
218  ON_Surface **t3,
219  ON_Surface **t4)
220 {
221  if (!(*t1)) {
222  ON_NurbsSurface *nt1 = ON_NurbsSurface::New();
223  (*t1)= (ON_Surface *)(nt1);
224  }
225  if (!(*t2)) {
226  ON_NurbsSurface *nt2 = ON_NurbsSurface::New();
227  (*t2)= (ON_Surface *)(nt2);
228  }
229  if (!(*t3)) {
230  ON_NurbsSurface *nt3 = ON_NurbsSurface::New();
231  (*t3)= (ON_Surface *)(nt3);
232  }
233  if (!(*t4)) {
234  ON_NurbsSurface *nt4 = ON_NurbsSurface::New();
235  (*t4)= (ON_Surface *)(nt4);
236  }
237 }
238 
240  const ON_Surface *srf,
241  ON_Interval *u_val,
242  ON_Interval *v_val,
243  ON_Surface **t1,
244  ON_Surface **t2,
245  ON_Surface **t3,
246  ON_Surface **t4,
247  ON_Surface **result
248  )
249 {
250  bool split = true;
251  ON_Surface **target;
252  int last_split = 0;
253  int t1_del, t2_del, t3_del;
254 
255  // Make sure we have intervals with non-zero lengths
256  if ((u_val->Length() <= ON_ZERO_TOLERANCE) || (v_val->Length() <= ON_ZERO_TOLERANCE))
257  return false;
258 
259  // If we have the original surface domain, just return true
260  if ((fabs(u_val->Min() - srf->Domain(0).m_t[0]) <= ON_ZERO_TOLERANCE) &&
261  (fabs(u_val->Max() - srf->Domain(0).m_t[1]) <= ON_ZERO_TOLERANCE) &&
262  (fabs(v_val->Min() - srf->Domain(1).m_t[0]) <= ON_ZERO_TOLERANCE) &&
263  (fabs(v_val->Max() - srf->Domain(1).m_t[1]) <= ON_ZERO_TOLERANCE)) {
264  (*result) = (ON_Surface *)srf;
265  return true;
266  }
267  if (fabs(u_val->Min() - srf->Domain(0).m_t[0]) > ON_ZERO_TOLERANCE) last_split = 1;
268  if (fabs(u_val->Max() - srf->Domain(0).m_t[1]) > ON_ZERO_TOLERANCE) last_split = 2;
269  if (fabs(v_val->Min() - srf->Domain(1).m_t[0]) > ON_ZERO_TOLERANCE) last_split = 3;
270  if (fabs(v_val->Max() - srf->Domain(1).m_t[1]) > ON_ZERO_TOLERANCE) last_split = 4;
271  t1_del = (*t1) ? (0) : (1);
272  t2_del = (*t2) ? (0) : (1);
273  t3_del = (*t3) ? (0) : (1);
274  ON_Surface *ssplit = (ON_Surface *)srf;
275  ON_Surface_Create_Scratch_Surfaces(t1, t2, t3, t4);
276  if (fabs(u_val->Min() - srf->Domain(0).m_t[0]) > ON_ZERO_TOLERANCE) {
277  if (last_split == 1) {target = t4;} else {target = t2;}
278  split = ssplit->Split(0, u_val->Min(), *t1, *target);
279  ssplit = *target;
280  }
281  if ((fabs(u_val->Max() - srf->Domain(0).m_t[1]) > ON_ZERO_TOLERANCE) && split) {
282  if (last_split == 2) {target = t4;} else {target = t1;}
283  split = ssplit->Split(0, u_val->Max(), *target, *t3);
284  ssplit = *target;
285  }
286  if ((fabs(v_val->Min() - srf->Domain(1).m_t[0]) > ON_ZERO_TOLERANCE) && split) {
287  if (last_split == 3) {target = t4;} else {target = t3;}
288  split = ssplit->Split(1, v_val->Min(), *t2, *target);
289  ssplit = *target;
290  }
291  if ((fabs(v_val->Max() - srf->Domain(1).m_t[1]) > ON_ZERO_TOLERANCE) && split) {
292  split = ssplit->Split(1, v_val->Max(), *t4, *t2);
293  }
294  (*result) = *t4;
295  if (t1_del) delete *t1;
296  if (t2_del) delete *t2;
297  if (t3_del) delete *t3;
298  (*result)->SetDomain(0,u_val->Min(), u_val->Max());
299  (*result)->SetDomain(1,v_val->Min(), v_val->Max());
300  return split;
301 }
302 
304  const ON_Surface *surf,
305  const ON_Interval& u,
306  const ON_Interval& v,
307  double upt,
308  double vpt,
309  ON_Surface **q0,
310  ON_Surface **q1,
311  ON_Surface **q2,
312  ON_Surface **q3)
313 {
314  bool split_success = true;
315  ON_Surface *north = NULL;
316  ON_Surface *south = NULL;
317 
318  // upt and vpt must be within their respective domains
319  if (!u.Includes(upt, true) || !v.Includes(vpt, true)) return false;
320 
321  // All four output surfaces should be NULL - the point of this function is to create them
322  if ((*q0) || (*q1) || (*q2) || (*q3)) {
323  bu_log("ON_Surface_Quad_Split was supplied non-NULL surfaces as output targets: q0: %p, q1: %p, q2: %p, q3: %p\n",
324  static_cast<void *>(*q0), static_cast<void *>(*q1), static_cast<void *>(*q2), static_cast<void *>(*q3));
325  return false;
326  }
327 
328  // First, get the north and south pieces
329  split_success = surf->Split(1, vpt, south, north);
330  if (!split_success || !south || !north) {
331  delete south;
332  delete north;
333  return false;
334  }
335 
336  // Split the south pieces to get q0 and q1
337  split_success = south->Split(0, upt, (*q0), (*q1));
338  if (!split_success || !(*q0) || !(*q1)) {
339  delete south;
340  delete north;
341  if (*q0) delete (*q0);
342  if (*q1) delete (*q1);
343  return false;
344  }
345 
346  // Split the north pieces to get q2 and q3
347  split_success = north->Split(0, upt, (*q3), (*q2));
348  if (!split_success || !(*q3) || !(*q2)) {
349  delete south;
350  delete north;
351  if (*q0) delete (*q0);
352  if (*q1) delete (*q1);
353  if (*q2) delete (*q2);
354  if (*q3) delete (*q3);
355  return false;
356  }
357 
358  delete south;
359  delete north;
360 
361  return true;
362 }
363 
364 // Local Variables:
365 // tab-width: 8
366 // mode: C++
367 // c-basic-offset: 4
368 // indent-tabs-mode: t
369 // c-file-style: "stroustrup"
370 // End:
371 // ex: shiftwidth=4 tabstop=8
Definition: db_flip.c:35
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
bool ON_Surface_IsFlat_U(ON_Plane *frames, double f_tol)
Perform flatness test of surface in U only.
Header file for the BRL-CAD common definitions.
bool ON_Surface_SubSurface(const ON_Surface *srf, ON_Interval *u_val, ON_Interval *v_val, ON_Surface **t1, ON_Surface **t2, ON_Surface **t3, ON_Surface **t4, ON_Surface **result)
Create a surface based on a subset of a parent surface.
void ON_Surface_Create_Scratch_Surfaces(ON_Surface **t1, ON_Surface **t2, ON_Surface **t3, ON_Surface **t4)
Create surfaces and store their pointers in the t* arguments.
double ON_Curve_Get_Horizontal_Tangent(const ON_Curve *curve, double min, double max, double zero_tol)
Search for a horizontal tangent on the curve between two curve parameters.
double ON_Curve_Get_Tangent(int direction, const ON_Curve *curve, double min, double max, double zero_tol)
double ON_Curve_Get_Vertical_Tangent(const ON_Curve *curve, double min, double max, double zero_tol)
Search for a vertical tangent on the curve between two curve parameters.
bool ON_Surface_IsFlat_V(ON_Plane *frames, double f_tol)
Perform flatness test of surface in V only.
bool ON_Surface_IsStraight(ON_Plane *frames, double s_tol)
Perform straightness test of surface.
int ON_Curve_Has_Tangent(const ON_Curve *curve, double ct_min, double ct_max, double t_tol)
Test whether a curve interval contains one or more horizontal or vertical tangents.
bool ON_Surface_IsFlat(ON_Plane *frames, double f_tol)
Perform flatness test of surface.
bool ON_Surface_Quad_Split(const ON_Surface *surf, const ON_Interval &u, const ON_Interval &v, double upt, double vpt, ON_Surface **q0, ON_Surface **q1, ON_Surface **q2, ON_Surface **q3)
Create four sub-surfaces from a parent surface.
bool ON_NearZero(double val, double epsilon)
Return truthfully whether a value is within a specified epsilon distance from zero.