BRL-CAD
shape_recognition_cylinder.cpp
Go to the documentation of this file.
1 #include "common.h"
2 
3 #include <set>
4 #include <map>
5 
6 #include "bu/log.h"
7 #include "bu/str.h"
8 #include "bu/malloc.h"
9 #include "shape_recognition.h"
10 
11 // Remove degenerate edge sets. A degenerate edge set is defined as two
12 // linear segments having the same two vertices. (To be sure, we should probably
13 // check curve directions in loops in some fashion...)
14 void
16  std::set<int> degenerate;
17  std::set<int>::iterator e_it;
18  for (e_it = edges->begin(); e_it != edges->end(); e_it++) {
19  if (degenerate.find(*e_it) == degenerate.end()) {
20  ON_BrepEdge& edge = data->brep->m_E[*e_it];
21  if (edge.EdgeCurveOf()->IsLinear()) {
22  for (int j = 0; j < data->edges_cnt; j++) {
23  int f_ind = data->edges[j];
24  ON_BrepEdge& edge2 = data->brep->m_E[f_ind];
25  if (edge2.EdgeCurveOf()->IsLinear()) {
26  if ((edge.Vertex(0)->Point() == edge2.Vertex(0)->Point() && edge.Vertex(1)->Point() == edge2.Vertex(1)->Point()) ||
27  (edge.Vertex(1)->Point() == edge2.Vertex(0)->Point() && edge.Vertex(0)->Point() == edge2.Vertex(1)->Point()))
28  {
29  degenerate.insert(*e_it);
30  degenerate.insert(f_ind);
31  break;
32  }
33  }
34  }
35  }
36  }
37  }
38  for (e_it = degenerate.begin(); e_it != degenerate.end(); e_it++) {
39  //std::cout << "erasing " << *e_it << "\n";
40  edges->erase(*e_it);
41  }
42 }
43 
44 
45 int
47 {
48  std::set<int>::iterator f_it;
49  std::set<int> planar_surfaces;
50  std::set<int> cylindrical_surfaces;
51  // First, check surfaces. If a surface is anything other than a plane or cylindrical,
52  // the verdict is no. If we don't have at least two planar surfaces and one
53  // cylindrical, the verdict is no.
54  for (int i = 0; i < data->faces_cnt; i++) {
55  int f_ind = data->faces[i];
56  ON_BrepFace *used_face = &(data->brep->m_F[f_ind]);
57  int surface_type = (int)GetSurfaceType(used_face->SurfaceOf(), NULL);
58  switch (surface_type) {
59  case SURFACE_PLANE:
60  planar_surfaces.insert(f_ind);
61  break;
62  case SURFACE_CYLINDER:
63  cylindrical_surfaces.insert(f_ind);
64  break;
65  default:
66  return 0;
67  break;
68  }
69  }
70  if (planar_surfaces.size() < 2) return 0;
71  if (cylindrical_surfaces.size() < 1) return 0;
72 
73  // Second, check if all cylindrical surfaces share the same axis.
74  ON_Cylinder cylinder;
75  ON_Surface *cs = data->brep->m_F[*cylindrical_surfaces.begin()].SurfaceOf()->Duplicate();
76  cs->IsCylinder(&cylinder);
77  delete cs;
78  for (f_it = cylindrical_surfaces.begin(); f_it != cylindrical_surfaces.end(); f_it++) {
79  ON_Cylinder f_cylinder;
80  ON_Surface *fcs = data->brep->m_F[(*f_it)].SurfaceOf()->Duplicate();
81  fcs->IsCylinder(&f_cylinder);
82  delete fcs;
83  if (f_cylinder.circle.Center().DistanceTo(cylinder.circle.Center()) > BREP_CYLINDRICAL_TOL) return 0;
84  }
85  // Third, see if all planes are coplanar with two and only two planes.
86  ON_Plane p1, p2;
87  int p2_set = 0;
88  data->brep->m_F[*planar_surfaces.begin()].SurfaceOf()->IsPlanar(&p1);
89  for (f_it = planar_surfaces.begin(); f_it != planar_surfaces.end(); f_it++) {
90  ON_Plane f_p;
91  data->brep->m_F[(*f_it)].SurfaceOf()->IsPlanar(&f_p);
92  if (!p2_set && f_p != p1) {
93  p2 = f_p;
94  continue;
95  };
96  if (f_p != p1 && f_p != p2) return 0;
97  }
98 
99  // Fourth, check that the two planes are parallel to each other.
100  if (p1.Normal().IsParallelTo(p2.Normal(), cyl_tol) == 0) {
101  std::cout << "p1 Normal: " << p1.Normal().x << "," << p1.Normal().y << "," << p1.Normal().z << "\n";
102  std::cout << "p2 Normal: " << p2.Normal().x << "," << p2.Normal().y << "," << p2.Normal().z << "\n";
103  return 0;
104  }
105 
106  // Fifth, remove degenerate edge sets.
107  std::set<int> active_edges;
108  array_to_set(&active_edges, data->edges, data->edges_cnt);
109  subbrep_remove_degenerate_edges(data, &active_edges);
110 
111  // Sixth, check for any remaining linear segments. For rpc primitives
112  // those are expected, but for a true cylinder the linear segments should
113  // all wash out in the degenerate pass
114  std::set<int>::iterator e_it;
115  for (e_it = active_edges.begin(); e_it != active_edges.end(); e_it++) {
116  ON_BrepEdge& edge = data->brep->m_E[*e_it];
117  if (edge.EdgeCurveOf()->IsLinear()) return 0;
118  }
119 
120  // Seventh, sort the curved edges into one of two circles. Again, in more
121  // general cases we might have other curves but a true cylinder should have
122  // all of its arcs on these two circles. We don't need to check for closed
123  // loops because we are assuming that in the input Brep; any curve except
124  // arc curves that survived the degeneracy test has already resulted in an
125  // earlier rejection.
126  std::set<int> arc_set_1, arc_set_2;
127  ON_Circle set1_c, set2_c;
128  int arc1_circle_set= 0;
129  int arc2_circle_set = 0;
130  for (e_it = active_edges.begin(); e_it != active_edges.end(); e_it++) {
131  ON_BrepEdge& edge = data->brep->m_E[*e_it];
132  ON_Arc arc;
133  if (edge.EdgeCurveOf()->IsArc(NULL, &arc, cyl_tol)) {
134  int assigned = 0;
135  ON_Circle circ(arc.StartPoint(), arc.MidPoint(), arc.EndPoint());
136  //std::cout << "circ " << circ.Center().x << " " << circ.Center().y << " " << circ.Center().z << "\n";
137  if (!arc1_circle_set) {
138  arc1_circle_set = 1;
139  set1_c = circ;
140  //std::cout << "center 1 " << set1_c.Center().x << " " << set1_c.Center().y << " " << set1_c.Center().z << "\n";
141  } else {
142  if (!arc2_circle_set) {
143  if (!(NEAR_ZERO(circ.Center().DistanceTo(set1_c.Center()), cyl_tol))){
144  arc2_circle_set = 1;
145  set2_c = circ;
146  //std::cout << "center 2 " << set2_c.Center().x << " " << set2_c.Center().y << " " << set2_c.Center().z << "\n";
147  }
148  }
149  }
150  if (NEAR_ZERO(circ.Center().DistanceTo(set1_c.Center()), cyl_tol)){
151  arc_set_1.insert(*e_it);
152  assigned = 1;
153  }
154  if (arc2_circle_set) {
155  if (NEAR_ZERO(circ.Center().DistanceTo(set2_c.Center()), cyl_tol)){
156  arc_set_2.insert(*e_it);
157  assigned = 1;
158  }
159  }
160  if (!assigned) {
161  std::cout << "found extra circle - no go\n";
162  return 0;
163  }
164  }
165  }
166 
167  data->type = CYLINDER;
168 
169  ON_3dVector hvect(set2_c.Center() - set1_c.Center());
170 
171  data->params->origin[0] = set1_c.Center().x;
172  data->params->origin[1] = set1_c.Center().y;
173  data->params->origin[2] = set1_c.Center().z;
174  data->params->hv[0] = hvect.x;
175  data->params->hv[1] = hvect.y;
176  data->params->hv[2] = hvect.z;
177  data->params->radius = set1_c.Radius();
178 
179  return 1;
180 }
181 
182 
183 int
184 cylindrical_loop_planar_vertices(ON_BrepFace *face, int loop_index)
185 {
186  std::set<int> verts;
187  ON_Brep *brep = face->Brep();
188  ON_BrepLoop *loop = &(brep->m_L[loop_index]);
189  for (int ti = 0; ti < loop->m_ti.Count(); ti++) {
190  ON_BrepTrim& trim = brep->m_T[loop->m_ti[ti]];
191  if (trim.m_ei != -1) {
192  ON_BrepEdge& edge = brep->m_E[trim.m_ei];
193  verts.insert(edge.Vertex(0)->m_vertex_index);
194  verts.insert(edge.Vertex(1)->m_vertex_index);
195  }
196  }
197  if (verts.size() == 3) {
198  //std::cout << "Three points - planar.\n";
199  return 1;
200  } else if (verts.size() >= 3) {
201  std::set<int>::iterator v_it = verts.begin();
202  ON_3dPoint p1 = brep->m_V[*v_it].Point();
203  v_it++;
204  ON_3dPoint p2 = brep->m_V[*v_it].Point();
205  v_it++;
206  ON_3dPoint p3 = brep->m_V[*v_it].Point();
207  ON_Plane test_plane(p1, p2, p3);
208  for (v_it = verts.begin(); v_it != verts.end(); v_it++) {
209  if (!NEAR_ZERO(test_plane.DistanceTo(brep->m_V[*v_it].Point()), BREP_PLANAR_TOL)) {
210  //std::cout << "vertex " << *v_it << " too far from plane, not planar: " << test_plane.DistanceTo(brep->m_V[*v_it].Point()) << "\n";
211  return 0;
212  }
213  }
214  //std::cout << verts.size() << " points, planar\n";
215  return 1;
216  } else {
217  //std::cout << "Closed single curve loop - planar only if surface is.";
218  return 0;
219  }
220  return 0;
221 }
222 
223 int
225 {
226  std::set<int> loops;
227  std::set<int>::iterator l_it;
228  ON_BrepFace *face = &(data->brep->m_F[face_index]);
229  array_to_set(&loops, data->loops, data->loops_cnt);
230  for(l_it = loops.begin(); l_it != loops.end(); l_it++) {
231  return cylindrical_loop_planar_vertices(face, face->m_li[*l_it]);
232  }
233  return 0;
234 }
235 
236 int
238 {
239  bu_log("process partial cylinder\n");
240  std::set<int> planar_surfaces;
241  std::set<int> cylindrical_surfaces;
242  for (int i = 0; i < data->faces_cnt; i++) {
243  int f_ind = data->faces[i];
244  ON_BrepFace *used_face = &(data->brep->m_F[f_ind]);
245  int surface_type = (int)GetSurfaceType(used_face->SurfaceOf(), NULL);
246  switch (surface_type) {
247  case SURFACE_PLANE:
248  planar_surfaces.insert(f_ind);
249  break;
250  case SURFACE_CYLINDER:
251  cylindrical_surfaces.insert(f_ind);
252  break;
253  default:
254  std::cout << "what???\n";
255  return 0;
256  break;
257  }
258  }
259 
260  // Characterize the planes of the non-linear edges. We need two planes - more
261  // than that indicate some sort of subtraction behavior. TODO - Eventually we should
262  // be able to characterize which edges form the subtraction shape candidate, but
263  // for now just bail unless we're dealing with the simple case.
264  //
265  // TODO - this test is adequate only for RCC shapes. Need to generalize
266  // this to check for both arcs and shared planes in non-arc curves to
267  // accomidate csg situations.
268  std::set<int> arc_set_1, arc_set_2;
269  ON_Circle set1_c, set2_c;
270  int arc1_circle_set= 0;
271  int arc2_circle_set = 0;
272 
273  for (int i = 0; i < data->edges_cnt; i++) {
274  int ei = data->edges[i];
275  ON_BrepEdge& edge = data->brep->m_E[ei];
276  if (!edge.EdgeCurveOf()->IsLinear()) {
277 
278  ON_Arc arc;
279  if (edge.EdgeCurveOf()->IsArc(NULL, &arc, cyl_tol)) {
280  int assigned = 0;
281  ON_Circle circ(arc.StartPoint(), arc.MidPoint(), arc.EndPoint());
282  //std::cout << "circ " << circ.Center().x << " " << circ.Center().y << " " << circ.Center().z << "\n";
283  if (!arc1_circle_set) {
284  arc1_circle_set = 1;
285  set1_c = circ;
286  //std::cout << "center 1 " << set1_c.Center().x << " " << set1_c.Center().y << " " << set1_c.Center().z << "\n";
287  } else {
288  if (!arc2_circle_set) {
289  if (!(NEAR_ZERO(circ.Center().DistanceTo(set1_c.Center()), cyl_tol))){
290  arc2_circle_set = 1;
291  set2_c = circ;
292  //std::cout << "center 2 " << set2_c.Center().x << " " << set2_c.Center().y << " " << set2_c.Center().z << "\n";
293  }
294  }
295  }
296  if (NEAR_ZERO(circ.Center().DistanceTo(set1_c.Center()), cyl_tol)){
297  arc_set_1.insert(ei);
298  assigned = 1;
299  }
300  if (arc2_circle_set) {
301  if (NEAR_ZERO(circ.Center().DistanceTo(set2_c.Center()), cyl_tol)){
302  arc_set_2.insert(ei);
303  assigned = 1;
304  }
305  }
306  if (!assigned) {
307  std::cout << "found extra circle - no go\n";
308  return 0;
309  }
310  }
311  }
312  }
313 
314 
315  // CSG representable cylinders may represent one or both of the
316  // following cases:
317  //
318  // a) non-parallel end caps - one or both capping planes are not
319  // perpendicular to the axis of the cylinder.
320  //
321  // b) partial cylindrical surface - some portion of the cylinderical
322  // surface is trimmed away.
323  //
324  // There are an infinite number of ways in which subsets of a cylinder
325  // may be removed by trimming curves - the plan is for complexities
326  // introduced into the outer loops to be reduced by recognizing the
327  // complex portions of those curves as influences of other shapes.
328  // Once recognized, the loops are simplified until we reach a shape
329  // that can be handled by the above cases.
330  //
331  // For example, let's say a cylindrical face has the following
332  // trim loop in its UV space:
333  //
334  // -------------------------
335  // | |
336  // | |
337  // | ************* |
338  // | * * |
339  // ------- -------
340  //
341  // The starred portion of the trimming curve is not representable
342  // in this CSG scheme, but if that portion of the curve is the
343  // result of a subtraction of another shape in the parent brep,
344  // then that portion of the curve can be treated as implicit in
345  // the subtraction of that other object. The complex lower trim
346  // curve set can then be replaced by a line between the two corner
347  // vertex points, which are not removed by the subtraction.
348  //
349  // Until such cases can be resolved, any curve complications of
350  // this sort are a conversion blocker. To make sure the supplied
351  // inputs are cases that can be handled, we collect all of the
352  // vertices in the data that are connected to one and only one
353  // non-linear edge in the set. Failure cases are:
354  //
355  // * More than four vertices that are mated with exactly one
356  // non-linear edge in the data set
357  // * Four vertices meeting previous criteria that are non-planar
358  // * Any vertex on a linear edge that is not coplanar with the
359  // plane described by the vertices meeting the above criteria
360  std::set<int> candidate_verts;
361  std::set<int> corner_verts; /* verts with one nonlinear edge */
362  std::set<int> linear_verts; /* verts with only linear edges */
363  std::set<int>::iterator v_it, e_it;
364  std::set<int> edges;
365  array_to_set(&edges, data->edges, data->edges_cnt);
366  // collect all candidate vertices
367  for (int i = 0; i < data->edges_cnt; i++) {
368  int ei = data->edges[i];
369  ON_BrepEdge& edge = data->brep->m_E[ei];
370  candidate_verts.insert(edge.Vertex(0)->m_vertex_index);
371  candidate_verts.insert(edge.Vertex(1)->m_vertex_index);
372  }
373  for (v_it = candidate_verts.begin(); v_it != candidate_verts.end(); v_it++) {
374  ON_BrepVertex& vert = data->brep->m_V[*v_it];
375  int curve_cnt = 0;
376  int line_cnt = 0;
377  for (int i = 0; i < vert.m_ei.Count(); i++) {
378  int ei = vert.m_ei[i];
379  ON_BrepEdge& edge = data->brep->m_E[ei];
380  if (edges.find(edge.m_edge_index) != edges.end()) {
381  if (edge.EdgeCurveOf()->IsLinear()) {
382  line_cnt++;
383  } else {
384  curve_cnt++;
385  }
386  }
387  }
388  if (curve_cnt == 1) {
389  corner_verts.insert(*v_it);
390  //std::cout << "found corner vert: " << *v_it << "\n";
391  }
392  if (line_cnt > 1 && curve_cnt == 0) {
393  linear_verts.insert(*v_it);
394  std::cout << "found linear vert: " << *v_it << "\n";
395  }
396  }
397 
398  // First, check corner count
399  if (corner_verts.size() > 4) {
400  std::cout << "more than 4 corners - complex\n";
401  return 0;
402  }
403 
404  // Second, create the candidate face plane. Verify coplanar status of points if we've got 4.
405  ON_Plane pcyl;
406  if (corner_verts.size() == 4) {
407  std::set<int>::iterator s_it = corner_verts.begin();
408  ON_3dPoint p1 = data->brep->m_V[*v_it].Point();
409  s_it++;
410  ON_3dPoint p2 = data->brep->m_V[*v_it].Point();
411  s_it++;
412  ON_3dPoint p3 = data->brep->m_V[*v_it].Point();
413  s_it++;
414  ON_3dPoint p4 = data->brep->m_V[*v_it].Point();
415  ON_Plane tmp_plane(p1, p2, p3);
416  if (tmp_plane.DistanceTo(p4) > BREP_PLANAR_TOL) {
417  std::cout << "planar tol fail\n";
418  return 0;
419  } else {
420  pcyl = tmp_plane;
421  }
422  } else {
423  // TODO - If we have less than four corner points and no additional curve planes, we
424  // must have a face subtraction that tapers to a point at the edge of the
425  // cylinder. Pull the linear edges from the two corner points to find the third point -
426  // this is a situation where a simpler arb (arb6?) is adequate to make the subtraction.
427  }
428 
429  // Third, if we had vertices with only linear edges, check to make sure they are in fact
430  // on the candidate plane for the face (if not, we've got something more complex going on).
431  if (linear_verts.size() > 0) {
432  std::set<int>::iterator s_it;
433  for (s_it = linear_verts.begin(); s_it != linear_verts.end(); s_it++) {
434  ON_3dPoint pnt = data->brep->m_V[*s_it].Point();
435  if (pcyl.DistanceTo(pnt) > BREP_PLANAR_TOL) {
436  std::cout << "stray verts fail\n";
437  return 0;
438  }
439  }
440  }
441 
442 #if 0
443  // Now, get the cylindrical surface properties.
444  ON_Cylinder cylinder;
445  ON_Surface *cs = data->brep->m_F[*cylindrical_surfaces.begin()].SurfaceOf()->Duplicate();
446  cs->IsCylinder(&cylinder, cyl_tol);
447  delete cs;
448 #endif
449 
450  // Check if the two circles are parallel to each other. If they are, and we have
451  // no corner points, then we have a complete cylinder
452  if (set1_c.Plane().Normal().IsParallelTo(set2_c.Plane().Normal(), cyl_tol) != 0) {
453  if (corner_verts.size() == 0) {
454  std::cout << "Full cylinder\n";
455  data->type = CYLINDER;
456 
457  ON_3dVector hvect(set2_c.Center() - set1_c.Center());
458 
459  data->params->origin[0] = set1_c.Center().x;
460  data->params->origin[1] = set1_c.Center().y;
461  data->params->origin[2] = set1_c.Center().z;
462  data->params->hv[0] = hvect.x;
463  data->params->hv[1] = hvect.y;
464  data->params->hv[2] = hvect.z;
465  data->params->radius = set1_c.Radius();
466 
467  return 1;
468  } else {
469  // We have parallel faces and corners - we need to subtract an arb
470  data->type = COMB;
471  std::cout << "Minus body arb\n";
472  return 1;
473  }
474  } else {
475  if (corner_verts.size() == 0) {
476  // We have non parallel faces and no corners - at least one and possible
477  // both end caps need subtracting, but no other subtractions are needed.
478  data->type = COMB;
479  std::cout << "Minus one or more end-cap arbs\n";
480  return 1;
481  } else {
482  // We have non parallel faces and corners - at least one and possible
483  // both end caps need subtracting, plus an arb to remove part of the
484  // cylinder body.
485  data->type = COMB;
486  std::cout << "Minus one or more end-cap arbs and body arb\n";
487  return 1;
488  }
489  }
490 
491 }
492 
493 // Local Variables:
494 // tab-width: 8
495 // mode: C++
496 // c-basic-offset: 4
497 // indent-tabs-mode: t
498 // c-file-style: "stroustrup"
499 // End:
500 // ex: shiftwidth=4 tabstop=8
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
Header file for the BRL-CAD common definitions.
#define BREP_CYLINDRICAL_TOL
COMPLEX data[64]
Definition: fftest.c:34
int cylindrical_planar_vertices(struct subbrep_object_data *data, int face_index)
surface_t GetSurfaceType(const ON_Surface *in_surface, struct filter_obj *obj)
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
void array_to_set(std::set< int > *set, int *array, int array_cnt)
int cylinder_csg(struct subbrep_object_data *data, fastf_t cyl_tol)
Definition: joint.h:84
void subbrep_remove_degenerate_edges(struct subbrep_object_data *data, std::set< int > *edges)
HIDDEN const double cs
Definition: sh_prj.c:617
#define BREP_PLANAR_TOL
int subbrep_is_cylinder(struct subbrep_object_data *data, fastf_t cyl_tol)
double fastf_t
Definition: defines.h:300
int cylindrical_loop_planar_vertices(ON_BrepFace *face, int loop_index)
csg_object_params * params