BRL-CAD
test_bot2nurbs.cpp
Go to the documentation of this file.
1 /* T E S T _ B O T 2 N U R B 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 test_bot2nurbs.cpp
21  *
22  * This file contains the logic that takes an input BoT, breaks that
23  * BoT down into patches, fits surfaces to those patches, builds
24  * trimming loops for the surfaces based on NURBS curves fitted to
25  * patch edge segments, and assembles the result into a closed BREP.
26  *
27  */
28 
29 /* TODO - Can't keep shifting the fitting plane for a triangle patch
30  * every time a triangle is added - too computationally expensive.
31  * Instead, need to make sure that the starting plane used for the
32  * NURBS fit is the same plane being used for the triangle overlap
33  * testing. Hopefully, that will preclude any "vertical or worse"
34  * surface fit issues while still allowing for a good surface fit,
35  * and avoid this business of triangles that overlap in one fit
36  * plane but not in another. Will probably need a change or two
37  * in the PCL surface fitting code.*/
38 
39 
40 #include "common.h"
41 
42 #include <map>
43 #include <set>
44 #include <queue>
45 #include <list>
46 #include <vector>
47 #include <iostream>
48 #include <fstream>
49 
50 #include "vmath.h"
51 #include "raytrace.h"
52 #include "wdb.h"
53 #include "plot3.h"
54 #include "opennurbs.h"
55 #include "opennurbs_fit.h"
56 #include <Eigen/SVD>
57 #include "brep.h"
58 
59 typedef std::pair<size_t, size_t> Edge;
60 typedef std::map< size_t, std::set<Edge> > VertToEdge;
61 typedef std::map< size_t, std::set<size_t> > VertToPatch;
62 typedef std::map< Edge, std::set<size_t> > EdgeToPatch;
63 typedef std::map< Edge, std::set<size_t> > EdgeToFace;
64 typedef std::map< size_t, size_t > FaceToPatch;
65 typedef std::map< size_t, size_t > PatchToPlane;
66 typedef std::set<Edge> EdgeList;
67 typedef std::set<size_t> FaceList;
68 
69 struct Manifold_Info {
72  struct rt_bot_internal *bot;
73 
74  std::map< size_t, std::set<size_t> > patches;
75  std::map< size_t, std::set<size_t> > patch_edges;
76  std::map< size_t, std::vector<size_t> > polycurves;
77 
83  std::map< size_t, size_t > face_to_plane;
85  ON_3dVectorArray face_normals;
86  ON_3dVectorArray vectors;
87  std::map< std::pair<size_t, size_t>, fastf_t > norm_results;
91  int patch_cnt;
92 };
93 
94 
95 /**********************************************************************************
96  * Utility functions used throughout the process of fitting
97  **********************************************************************************/
98 
99 // Make an edge using consistent vertex ordering
100 static Edge
101 mk_edge(size_t pt_A, size_t pt_B)
102 {
103  if (pt_A <= pt_B) {
104  return std::make_pair(pt_A, pt_B);
105  } else {
106  return std::make_pair(pt_B, pt_A);
107  }
108 }
109 
110 
111 // Calculate area of a face
112 static double
113 face_area(struct rt_bot_internal *bot, size_t face_num)
114 {
115  point_t ptA, ptB, ptC;
116  double a, b, c, p;
117  double area;
118  VMOVE(ptA, &bot->vertices[bot->faces[face_num*3+0]*3]);
119  VMOVE(ptB, &bot->vertices[bot->faces[face_num*3+1]*3]);
120  VMOVE(ptC, &bot->vertices[bot->faces[face_num*3+2]*3]);
121  a = DIST_PT_PT(ptA, ptB);
122  b = DIST_PT_PT(ptB, ptC);
123  c = DIST_PT_PT(ptC, ptA);
124  p = (a + b + c)/2;
125  area = sqrt(p*(p-a)*(p-b)*(p-c));
126  return area;
127 }
128 
129 
130 // Given a patch consisting of a set of faces find the
131 // set of edges that forms the outer edge of the patch
132 static void
133 find_edge_segments(std::set<size_t> *faces, EdgeList *patch_edges, struct Manifold_Info *info)
134 {
135  size_t pt_A, pt_B, pt_C;
136  std::set<size_t>::iterator it;
137  std::map<std::pair<size_t, size_t>, size_t> edge_face_cnt;
138  std::map<std::pair<size_t, size_t>, size_t>::iterator efc_it;
139  patch_edges->clear();
140  for (it = faces->begin(); it != faces->end(); it++) {
141  pt_A = info->bot->faces[(*it)*3+0]*3;
142  pt_B = info->bot->faces[(*it)*3+1]*3;
143  pt_C = info->bot->faces[(*it)*3+2]*3;
144  edge_face_cnt[mk_edge(pt_A, pt_B)] += 1;
145  edge_face_cnt[mk_edge(pt_B, pt_C)] += 1;
146  edge_face_cnt[mk_edge(pt_C, pt_A)] += 1;
147  }
148 
149  for (efc_it = edge_face_cnt.begin(); efc_it!=edge_face_cnt.end(); efc_it++) {
150  if ((*efc_it).second == 1) {
151  info->vert_to_edge[(*efc_it).first.first].insert((*efc_it).first);
152  info->vert_to_edge[(*efc_it).first.second].insert((*efc_it).first);
153  patch_edges->insert((*efc_it).first);
154  }
155  }
156 }
157 
158 
159 // Given a BoT face, return the three edges associated with that face
160 static void
161 get_face_edges(struct rt_bot_internal *bot, size_t face_num, std::set<Edge> *face_edges)
162 {
163  size_t pt_A, pt_B, pt_C;
164  pt_A = bot->faces[face_num*3+0]*3;
165  pt_B = bot->faces[face_num*3+1]*3;
166  pt_C = bot->faces[face_num*3+2]*3;
167  face_edges->insert(mk_edge(pt_A, pt_B));
168  face_edges->insert(mk_edge(pt_B, pt_C));
169  face_edges->insert(mk_edge(pt_C, pt_A));
170 }
171 
172 
173 // Given a BoT face, return the three other faces that share an edge with that face
174 static void
175 get_connected_faces(struct rt_bot_internal *bot, size_t face_num, EdgeToFace *edge_to_face, std::set<size_t> *faces)
176 {
177  std::set<Edge> face_edges;
178  std::set<Edge>::iterator face_edges_it;
179  get_face_edges(bot, face_num, &face_edges);
180  for (face_edges_it = face_edges.begin(); face_edges_it != face_edges.end(); face_edges_it++) {
181  std::set<size_t> faces_from_edge = (*edge_to_face)[(*face_edges_it)];
182  std::set<size_t>::iterator ffe_it;
183  for (ffe_it = faces_from_edge.begin(); ffe_it != faces_from_edge.end() ; ffe_it++) {
184  if ((*ffe_it) != face_num) {
185  faces->insert((*ffe_it));
186  }
187  }
188  }
189 }
190 
191 
192 // Use SVD algorithm from Soderkvist to fit a plane to vertex points
193 // http://www.math.ltu.se/~jove/courses/mam208/svd.pdf
194 static void
195 fit_plane(size_t UNUSED(patch_id), std::set<size_t> *faces, struct Manifold_Info *info, ON_Plane *plane)
196 {
197  if (faces->size() > 0) {
198  ON_3dPoint center(0.0, 0.0, 0.0);
199  std::set<size_t> verts;
200  std::set<size_t>::iterator f_it, v_it;
201  for (f_it = faces->begin(); f_it != faces->end(); f_it++) {
202  verts.insert(info->bot->faces[(*f_it)*3+0]*3);
203  verts.insert(info->bot->faces[(*f_it)*3+1]*3);
204  verts.insert(info->bot->faces[(*f_it)*3+2]*3);
205  }
206  point_t pt;
207  for (v_it = verts.begin(); v_it != verts.end(); v_it++) {
208  VMOVE(pt, &info->bot->vertices[*v_it]);
209  center.x += pt[0]/verts.size();
210  center.y += pt[1]/verts.size();
211  center.z += pt[2]/verts.size();
212  }
213  Eigen::MatrixXd A(3, verts.size());
214  int vert_cnt = 0;
215  for (v_it = verts.begin(); v_it != verts.end(); v_it++) {
216  VMOVE(pt, &info->bot->vertices[*v_it]);
217  A(0,vert_cnt) = pt[0] - center.x;
218  A(1,vert_cnt) = pt[1] - center.y;
219  A(2,vert_cnt) = pt[2] - center.z;
220  vert_cnt++;
221  }
222 
223  Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU);
224 
225  // 4. Normal is in column 3 of U matrix
226  ON_3dVector normal(svd.matrixU()(0,2), svd.matrixU()(1,2), svd.matrixU()(2,2));
227 
228  // 5. Construct plane
229  ON_Plane new_plane(center, normal);
230  (*plane) = new_plane;
231 /* Uncomment to plot fitting planes */
232 #if 0
233  struct bu_vls name;
234  bu_vls_init(&name);
235  bu_vls_printf(&name, "fit_plane_%d.pl", patch_id);
236  FILE* plot_file = fopen(bu_vls_addr(&name), "w");
237  int r = int(256*drand48() + 1.0);
238  int g = int(256*drand48() + 1.0);
239  int b = int(256*drand48() + 1.0);
240  point_t pc;
241  for (f_it = faces->begin(); f_it != faces->end(); f_it++) {
242  point_t p1, p2, p3;
243  VMOVE(p1, &info->bot->vertices[info->bot->faces[(*f_it)*3+0]*3]);
244  VMOVE(p2, &info->bot->vertices[info->bot->faces[(*f_it)*3+1]*3]);
245  VMOVE(p3, &info->bot->vertices[info->bot->faces[(*f_it)*3+2]*3]);
246  VSUB2(p1, p1, pc);
247  VSUB2(p2, p2, pc);
248  VSUB2(p3, p3, pc);
249  pdv_3move(plot_file, p1);
250  pdv_3cont(plot_file, p2);
251  pdv_3move(plot_file, p1);
252  pdv_3cont(plot_file, p3);
253  pdv_3move(plot_file, p2);
254  pdv_3cont(plot_file, p3);
255  }
256  pl_color(plot_file, 255, 255, 255);
257  point_t pn;
258  VSET(pc, center.x,center.y,center.z);
259  pdv_3move(plot_file, pc);
260  VSET(pn, normal.x,normal.y,normal.z);
261  VSCALE(pn, pn, 200);
262  VADD2(pn,pn,pc);
263  pdv_3cont(plot_file, pn);
264 
265  fclose(plot_file);
266 #endif
267  }
268 }
269 
270 
271 // Check if a mesh is planar
272 static int
273 planar_patch_test(size_t patch_id, std::set<size_t> *faces, struct Manifold_Info *info, fastf_t tol)
274 {
275  if (faces->size() == 1) return 1;
276  if (faces->size() > 0) {
277  std::set<size_t> verts;
278  std::set<size_t>::iterator f_it, v_it;
279  fastf_t max_dist = 0.0;
280  ON_Plane best_fit_plane;
281  fit_plane(patch_id, faces, info, &best_fit_plane);
282  for (f_it = faces->begin(); f_it != faces->end(); f_it++) {
283  verts.insert(info->bot->faces[(*f_it)*3+0]*3);
284  verts.insert(info->bot->faces[(*f_it)*3+1]*3);
285  verts.insert(info->bot->faces[(*f_it)*3+2]*3);
286  }
287  for (v_it = verts.begin(); v_it != verts.end(); v_it++) {
288  point_t pt;
289  VMOVE(pt, &info->bot->vertices[*v_it]);
290  fastf_t curr_dist = best_fit_plane.DistanceTo(ON_3dPoint(pt[0], pt[1], pt[2]));
291  if (curr_dist > max_dist) max_dist = curr_dist;
292  }
293  if (max_dist < tol) return 1;
294  return 0;
295  } else {
296  return 0;
297  }
298 }
299 
300 
301 /**********************************************************************************
302  * Debugging code for plotting structures commonly generated during fitting
303  **********************************************************************************/
304 
305 // To plot specific groups of curves in MGED, do the following:
306 // set glob_compat_mode 0
307 // set pl_list [glob 27_curve*.pl]
308 // foreach plfile $pl_list {overlay $plfile}
309 static void
310 plot_curve(struct rt_bot_internal *bot, std::vector<size_t> *pts, int r, int g, int b, FILE *c_plot)
311 {
312  if (c_plot == NULL) {
313  FILE* plot = fopen("plot_curve.pl", "w");
314  pl_color(plot, int(256*drand48() + 1.0), int(256*drand48() + 1.0), int(256*drand48() + 1.0));
315  std::vector<size_t>::iterator v_it;
316  for (v_it = pts->begin(); v_it != pts->end()-1; v_it++) {
317  pdv_3move(plot, &bot->vertices[(*v_it)]);
318  pdv_3cont(plot, &bot->vertices[(*(v_it+1))]);
319  }
320  fclose(plot);
321  } else {
322  pl_color(c_plot, r, g, b);
323  std::vector<size_t>::iterator v_it;
324  for (v_it = pts->begin(); v_it != pts->end()-1; v_it++) {
325  pdv_3move(c_plot, &bot->vertices[(*v_it)]);
326  pdv_3cont(c_plot, &bot->vertices[(*(v_it+1))]);
327  }
328  }
329 }
330 
331 
332 static void
333 plot_face(point_t p1, point_t p2, point_t p3, int r, int g, int b, FILE *c_plot)
334 {
335  pl_color(c_plot, r, g, b);
336  pdv_3move(c_plot, p1);
337  pdv_3cont(c_plot, p2);
338  pdv_3move(c_plot, p1);
339  pdv_3cont(c_plot, p3);
340  pdv_3move(c_plot, p2);
341  pdv_3cont(c_plot, p3);
342 }
343 
344 
345 static void
346 plot_faces(std::map< size_t, std::set<size_t> > *patches, struct Manifold_Info *info, const char *filename)
347 {
348  std::map< size_t, std::set<size_t> >::iterator p_it;
349  FILE* plot_file = fopen(filename, "w");
350  for (p_it = patches->begin(); p_it != patches->end(); p_it++) {
351  int r = int(256*drand48() + 1.0);
352  int g = int(256*drand48() + 1.0);
353  int b = int(256*drand48() + 1.0);
354  std::set<size_t> *faces = &((*p_it).second);
355  std::set<size_t>::iterator f_it;
356  for (f_it = faces->begin(); f_it != faces->end(); f_it++) {
357  plot_face(&info->bot->vertices[info->bot->faces[(*f_it)*3+0]*3], &info->bot->vertices[info->bot->faces[(*f_it)*3+1]*3], &info->bot->vertices[info->bot->faces[(*f_it)*3+2]*3], r, g ,b, plot_file);
358  }
359  }
360  fclose(plot_file);
361 }
362 
363 
364 static void
365 plot_patch_borders(std::map< size_t, std::set<size_t> > *patches, struct Manifold_Info *info, const char *filename)
366 {
367  std::map< size_t, std::set<size_t> >::iterator p_it;
368  FILE* edge_plot = fopen(filename, "w");
369  pl_color(edge_plot, 255, 255, 255);
370  for (p_it = patches->begin(); p_it != patches->end(); p_it++) {
371  EdgeList patch_edges;
372  EdgeList::iterator e_it;
373  find_edge_segments(&((*p_it).second), &patch_edges, info);
374  for (e_it = patch_edges.begin(); e_it != patch_edges.end(); e_it++) {
375  pdv_3move(edge_plot, &info->bot->vertices[(*e_it).first]);
376  pdv_3cont(edge_plot, &info->bot->vertices[(*e_it).second]);
377  }
378  }
379  fclose(edge_plot);
380 }
381 
382 
383 static void
384 plot_ncurve(ON_Curve &curve, FILE *c_plot)
385 {
386  double pt1[3], pt2[3];
387  ON_2dPoint from, to;
388  int plotres = 100;
389  if (curve.IsLinear()) {
390  int knotcnt = curve.SpanCount();
391  double *knots = new double[knotcnt + 1];
392 
393  curve.GetSpanVector(knots);
394  for (int i = 1; i <= knotcnt; i++) {
395  ON_3dPoint p = curve.PointAt(knots[i - 1]);
396  VMOVE(pt1, p);
397  p = curve.PointAt(knots[i]);
398  VMOVE(pt2, p);
399  pdv_3move(c_plot, pt1);
400  pdv_3cont(c_plot, pt2);
401  }
402 
403  } else {
404  ON_Interval dom = curve.Domain();
405  // XXX todo: dynamically sample the curve
406  for (int i = 1; i <= plotres; i++) {
407  ON_3dPoint p = curve.PointAt(dom.ParameterAt((double)(i - 1)
408  / (double)plotres));
409  VMOVE(pt1, p);
410  p = curve.PointAt(dom.ParameterAt((double) i / (double)plotres));
411  VMOVE(pt2, p);
412  pdv_3move(c_plot, pt1);
413  pdv_3cont(c_plot, pt2);
414  }
415  }
416 }
417 
418 
419 // plot loop
420 static void
421 plot_loop(std::vector<size_t> *edges, struct Manifold_Info *info, FILE *l_plot)
422 {
423  std::vector<size_t>::iterator v_it;
424  int r = int(256*drand48() + 1.0);
425  int g = int(256*drand48() + 1.0);
426  int b = int(256*drand48() + 1.0);
427  for (v_it = edges->begin(); v_it != edges->end(); v_it++) {
428  pl_color(l_plot, r, g, b);
429  ON_BrepEdge& edge = info->brep->m_E[(int)(*v_it)];
430  const ON_Curve* edge_curve = edge.EdgeCurveOf();
431  ON_Interval dom = edge_curve->Domain();
432  // XXX todo: dynamically sample the curve
433  int plotres = 50;
434  for (int i = 1; i <= plotres; i++) {
435  double pt1[3], pt2[3];
436  ON_3dPoint p = edge_curve->PointAt(dom.ParameterAt((double)(i - 1)
437  / (double)plotres));
438  VMOVE(pt1, p);
439  p = edge_curve->PointAt(dom.ParameterAt((double) i / (double)plotres));
440  VMOVE(pt2, p);
441  pdv_3move(l_plot, pt1);
442  pdv_3cont(l_plot, pt2);
443  }
444  }
445 }
446 
447 
448 /**********************************************************************************
449  * Code for partitioning a mesh into patches suitable for NURBS surface fitting
450  **********************************************************************************/
451 
452 // Clean out empty patches
453 static void
454 remove_empty_patches(struct Manifold_Info *info)
455 {
456  for (int i = 0; i < (int)info->patches.size(); i++) {
457  if (info->patches[i].size() == 0) {
458  info->patches.erase(i);
459  }
460  }
461 }
462 
463 
464 static void
465 sync_structure_maps(struct Manifold_Info *info)
466 {
467  std::map< size_t, std::set<size_t> >::iterator p_it;
468  info->edge_to_patch.clear();
469  info->vert_to_patch.clear();
470  info->face_to_patch.clear();
471  for (p_it = info->patches.begin(); p_it != info->patches.end(); p_it++) {
472  std::set<size_t>::iterator f_it;
473  std::set<size_t> *faces = &((*p_it).second);
474  for (f_it = faces->begin(); f_it != faces->end(); f_it++) {
475  info->face_to_patch[(*f_it)] = (*p_it).first;
476  std::set<Edge> face_edges;
477  std::set<Edge>::iterator face_edges_it;
478  get_face_edges(info->bot, (*f_it), &face_edges);
479  for (face_edges_it = face_edges.begin(); face_edges_it != face_edges.end(); face_edges_it++) {
480  info->edge_to_patch[(*face_edges_it)].insert((*p_it).first);
481  info->vert_to_patch[(*face_edges_it).first].insert((*p_it).first);
482  info->vert_to_patch[(*face_edges_it).second].insert((*p_it).first);
483  }
484  }
485  info->patch_to_plane[(*p_it).first] = info->face_to_plane[*(faces->begin())];
486  }
487 }
488 
489 
490 // For a given face, determine which patch shares the majority of its edges. If no one patch
491 // has a majority, return -1
492 static int
493 find_major_patch(size_t face_num, struct Manifold_Info *info, size_t curr_patch, size_t patch_count)
494 {
495  std::map<size_t, size_t> patch_cnt;
496  std::map<size_t, size_t>::iterator patch_cnt_itr;
497  std::set<size_t> connected_faces;
498  std::set<size_t>::iterator cf_it;
499  size_t max_patch = 0;
500  size_t max_patch_edge_cnt = 0;
501  get_connected_faces(info->bot, face_num, &(info->edge_to_face), &connected_faces);
502  for (cf_it = connected_faces.begin(); cf_it != connected_faces.end() ; cf_it++) {
503  patch_cnt[info->face_to_patch[*cf_it]] += 1;
504  }
505  for (patch_cnt_itr = patch_cnt.begin(); patch_cnt_itr != patch_cnt.end() ; patch_cnt_itr++) {
506  if ((*patch_cnt_itr).second > max_patch_edge_cnt) {
507  max_patch_edge_cnt = (*patch_cnt_itr).second;
508  max_patch = (*patch_cnt_itr).first;
509  }
510  }
511  if (max_patch_edge_cnt == 1) {
512  if (patch_count > 3) {
513  return -1;
514  } else {
515  // A patch with three or fewer triangles is considered a small patch - if one of the shared
516  // faces is a viable candidate, return it.
517  int candidate_patch = -1;
518  double candidate_patch_norm = 0;
519  for (cf_it = connected_faces.begin(); cf_it != connected_faces.end() ; cf_it++) {
520  if (info->face_to_patch[*cf_it] != curr_patch) {
521  double current_norm = info->norm_results[std::make_pair(face_num, info->patch_to_plane[info->face_to_patch[(*cf_it)]])];
522  if (current_norm > candidate_patch_norm) {
523  candidate_patch = info->face_to_patch[(*cf_it)];
524  candidate_patch_norm = current_norm;
525  }
526  }
527  }
528  return candidate_patch;
529  }
530  } else {
531  return (int)max_patch;
532  }
533 }
534 
535 
536 // Given a list of edges and the "current" patch, find triangles that share more edges with
537 // a different patch. If the corresponding normals allow, move the triangles in question to
538 // the patch with which they share the most edges. Return the number of triangles shifted.
539 // This won't fully "smooth" an edge in a given pass, but an iterative approach should converge
540 // to edges with triangles in their correct "major" patches
541 static size_t
542 shift_edge_triangles(std::map< size_t, std::set<size_t> > *patches, size_t curr_patch, struct Manifold_Info *info)
543 {
544  std::set<size_t> patch_edgefaces;
545  std::map<size_t, size_t> faces_to_shift;
546  std::map<size_t, size_t>::iterator f_it;
547  EdgeList edges;
548  EdgeList::iterator e_it;
549  std::set<size_t>::iterator sizet_it;
550 
551  // 0. Build edge set;
552  find_edge_segments(&((*patches)[curr_patch]), &edges, info);
553 
554  // 1. build triangle set of edge triangles in patch.
555  for (e_it = edges.begin(); e_it != edges.end(); e_it++) {
556  std::set<size_t> faces_from_edge = info->edge_to_face[(*e_it)];
557  for (sizet_it = faces_from_edge.begin(); sizet_it != faces_from_edge.end() ; sizet_it++) {
558  if (info->face_to_patch[(*sizet_it)] == curr_patch) {
559  patch_edgefaces.insert((*sizet_it));
560  }
561  }
562 
563  }
564  // 2. build triangle set of edge triangles with major patch other than current patch
565  // and a normal that permits movement to the major patch
566  for (sizet_it = patch_edgefaces.begin(); sizet_it != patch_edgefaces.end() ; sizet_it++) {
567  int major_patch = find_major_patch((*sizet_it), info, curr_patch, (*patches)[curr_patch].size());
568  if (major_patch != (int)curr_patch && major_patch != -1) {
569  if (info->norm_results[std::make_pair((*sizet_it), info->patch_to_plane[major_patch])] >= 0) {
570  faces_to_shift[(*sizet_it)] = major_patch;
571  }
572  }
573  }
574 
575  // 3. move triangles to their major patch.
576  for (f_it = faces_to_shift.begin(); f_it != faces_to_shift.end(); f_it++) {
577  (*patches)[info->face_to_patch[(*f_it).first]].erase((*f_it).first);
578  (*patches)[(*f_it).second].insert((*f_it).first);
579  info->face_to_patch[(*f_it).first] = (*f_it).second;
580  }
581 
582  return faces_to_shift.size();
583 }
584 
585 
586 static void
587 construct_patches(std::set<size_t> *faces, std::map< size_t, std::set<size_t> > *patches, struct Manifold_Info *info)
588 {
589  if (faces->size() > 0) {
590  while (faces->size() > 0) {
591  info->patch_cnt++;
592  size_t new_patch = info->patch_cnt;
593  std::queue<size_t> face_queue;
594  face_queue.push((*faces->begin()));
595  faces->erase(face_queue.front());
596  size_t start_face_num = face_queue.front();
597  info->patch_to_plane[new_patch] = info->face_to_plane[face_queue.front()];
598  while (!face_queue.empty()) {
599  size_t face_num = face_queue.front();
600  face_queue.pop();
601  (*patches)[new_patch].insert(face_num);
602  info->face_to_patch[face_num] = new_patch;
603  std::set<size_t> connected_faces;
604  std::set<size_t>::iterator cf_it;
605  get_connected_faces(info->bot, face_num, &(info->edge_to_face), &connected_faces);
606  for (cf_it = connected_faces.begin(); cf_it != connected_faces.end() ; cf_it++) {
607  if (faces->find((*cf_it)) != faces->end()) {
608  if (ON_DotProduct( *(*info).face_normals.At((int)(start_face_num)),*(*info).face_normals.At((int)(*cf_it)) ) > 0.5) {
609  face_queue.push((*cf_it));
610  faces->erase((*cf_it));
611  }
612  }
613  }
614  }
615  }
616  }
617 }
618 
619 
620 static void
621 split_overlapping_patch(size_t face1, size_t face2, size_t orig_patch, std::map< size_t, std::set<size_t> > *patches, struct Manifold_Info *info)
622 {
623  std::set<size_t> *faces = &((*patches)[orig_patch]);
624  std::queue<size_t> face_queue_1;
625  std::queue<size_t> face_queue_2;
626  face_queue_1.push(*faces->find(face1));
627  faces->erase(face_queue_1.front());
628  info->patch_cnt++;
629  size_t new_patch_1 = info->patch_cnt;
630  size_t start_face_num_1 = face_queue_1.front();
631  info->patch_to_plane[new_patch_1] = info->face_to_plane[face_queue_1.front()];
632  face_queue_2.push(*faces->find(face2));
633  faces->erase(face_queue_2.front());
634  info->patch_cnt++;
635  size_t new_patch_2 = info->patch_cnt;
636  size_t start_face_num_2 = face_queue_2.front();
637  info->patch_to_plane[new_patch_2] = info->face_to_plane[face_queue_2.front()];
638  while (!face_queue_1.empty() || !face_queue_2.empty()) {
639  if (!face_queue_1.empty()) {
640  size_t face_num_1 = face_queue_1.front();
641  face_queue_1.pop();
642  (*patches)[new_patch_1].insert(face_num_1);
643  info->face_to_patch[face_num_1] = new_patch_1;
644  std::set<size_t> connected_faces;
645  std::set<size_t>::iterator cf_it;
646  get_connected_faces(info->bot, face_num_1, &(info->edge_to_face), &connected_faces);
647  for (cf_it = connected_faces.begin(); cf_it != connected_faces.end() ; cf_it++) {
648  if (faces->find((*cf_it)) != faces->end() && (*patches)[new_patch_2].find(*cf_it) == (*patches)[new_patch_2].end()) {
649  if (ON_DotProduct( *(*info).face_normals.At((int)(start_face_num_1)),*(*info).face_normals.At((int)(*cf_it)) ) > 0.5) {
650  face_queue_1.push((*cf_it));
651  faces->erase((*cf_it));
652  }
653  }
654  }
655  }
656  if (!face_queue_2.empty()) {
657  size_t face_num_2 = face_queue_2.front();
658  face_queue_2.pop();
659  (*patches)[new_patch_2].insert(face_num_2);
660  info->face_to_patch[face_num_2] = new_patch_2;
661  std::set<size_t> connected_faces;
662  std::set<size_t>::iterator cf_it;
663  get_connected_faces(info->bot, face_num_2, &(info->edge_to_face), &connected_faces);
664  for (cf_it = connected_faces.begin(); cf_it != connected_faces.end(); cf_it++) {
665  if (faces->find((*cf_it)) != faces->end() && (*patches)[new_patch_1].find(*cf_it) == (*patches)[new_patch_1].end()) {
666  if (ON_DotProduct( *(*info).face_normals.At(((int)start_face_num_2)),*(*info).face_normals.At((int)(*cf_it)) ) > 0.5) {
667  face_queue_2.push((*cf_it));
668  faces->erase((*cf_it));
669  }
670  }
671  }
672  }
673  }
674  if (faces->size() > 0) construct_patches(faces, patches, info);
675 }
676 
677 
678 // Given a patch, find triangles that overlap when projected into the patch domain, remove them
679 // from the current patch and set them up in their own patch. Returns the number of faces moved
680 // from the current patch to their own patch.
681 //
682 // TODO - saw at least one case where triangle removal resulted one patch breaking into two
683 // topological patches. Need to handle this, otherwise Bad Things will happen when trying
684 // to generate BREPs.
685 // Also - need to check how the nurbs fitting code is calculating their initial plane, and use
686 // that for this test rather than the original faces - it is projection into *THAT* plane that
687 // must be clean, and a projection clean in the "standard" direction may not be clean there.
688 static size_t
689 overlapping_edge_triangles(std::map< size_t, std::set<size_t> > *patches, struct Manifold_Info *info)
690 {
691  size_t total_overlapping = 0;
692  std::map< size_t, std::set<size_t> >::iterator p_it;
693  for (p_it = patches->begin(); p_it != patches->end(); p_it++) {
694  if ((*p_it).second.size() > 0) {
695  ON_Plane plane;
696  fit_plane((*p_it).first, &((*p_it).second), info, &plane);
697  ON_Xform xf;
698  xf.PlanarProjection(plane);
699  size_t overlap_cnt = 1;
700  size_t patch_overlapping = 0;
701  while (overlap_cnt != 0) {
702  std::set<size_t> edge_triangles;
703  EdgeList edges;
704  find_edge_segments(&((*p_it).second), &edges, info);
705  EdgeList::iterator e_it;
706  std::set<size_t>::iterator sizet_it;
707  // 1. build triangle set of edge triangles in patch.
708  for (e_it = edges.begin(); e_it != edges.end(); e_it++) {
709  std::set<size_t> faces_from_edge = info->edge_to_face[(*e_it)];
710  for (sizet_it = faces_from_edge.begin(); sizet_it != faces_from_edge.end() ; sizet_it++) {
711  if (info->face_to_patch[(*sizet_it)] == (*p_it).first) {
712  edge_triangles.insert((*sizet_it));
713  }
714  }
715  }
716  // 2. Find an overlapping triangle, remove it to its own patch, fix the edges, continue
717  overlap_cnt = 0;
718  for (sizet_it = edge_triangles.begin(); sizet_it != edge_triangles.end() ; sizet_it++) {
719  point_t v0, v1, v2;
720  ON_3dPoint onpt_v0(&info->bot->vertices[info->bot->faces[(*sizet_it)*3+0]*3]);
721  ON_3dPoint onpt_v1(&info->bot->vertices[info->bot->faces[(*sizet_it)*3+1]*3]);
722  ON_3dPoint onpt_v2(&info->bot->vertices[info->bot->faces[(*sizet_it)*3+2]*3]);
723  onpt_v0.Transform(xf);
724  onpt_v1.Transform(xf);
725  onpt_v2.Transform(xf);
726  VSET(v0, onpt_v0.x, onpt_v0.y, onpt_v0.z);
727  VSET(v1, onpt_v1.x, onpt_v1.y, onpt_v1.z);
728  VSET(v2, onpt_v2.x, onpt_v2.y, onpt_v2.z);
729  edge_triangles.erase(*sizet_it);
730  std::set<size_t>::iterator ef_it2 = edge_triangles.begin();
731  for (ef_it2 = edge_triangles.begin(); ef_it2!=edge_triangles.end(); ef_it2++) {
732  if ((*sizet_it) != (*ef_it2)) {
733  point_t u0, u1, u2;
734  ON_3dPoint onpt_u0(&info->bot->vertices[info->bot->faces[(*ef_it2)*3+0]*3]);
735  ON_3dPoint onpt_u1(&info->bot->vertices[info->bot->faces[(*ef_it2)*3+1]*3]);
736  ON_3dPoint onpt_u2(&info->bot->vertices[info->bot->faces[(*ef_it2)*3+2]*3]);
737  onpt_u0.Transform(xf);
738  onpt_u1.Transform(xf);
739  onpt_u2.Transform(xf);
740  VSET(u0, onpt_u0.x, onpt_u0.y, onpt_u0.z);
741  VSET(u1, onpt_u1.x, onpt_u1.y, onpt_u1.z);
742  VSET(u2, onpt_u2.x, onpt_u2.y, onpt_u2.z);
743  int overlap = bn_tri_tri_isect_coplanar(v0, v1, v2, u0, u1, u2, 1);
744  if (overlap) {
745  std::cout << "Patch " << (*p_it).first << " overlap: (" << *sizet_it << "," << *ef_it2 << ")\n";
746  split_overlapping_patch(*sizet_it, *ef_it2, (*p_it).first, patches, info);
747 #if 0
748  info->patch_cnt++;
749  size_t new_patch = info->patch_cnt;
750  (*patches)[new_patch].insert(*sizet_it);
751  (*patches)[(*p_it).first].erase(*sizet_it);
752  info->patch_to_plane[new_patch] = info->face_to_plane[*sizet_it];
753  info->face_to_patch[*sizet_it] = new_patch;
754  find_edge_segments(&((*p_it).second), &edges, info);
755 #endif
756  overlap_cnt++;
757  patch_overlapping++;
758  total_overlapping++;
759  }
760  }
761  if (overlap_cnt > 0) break;
762  }
763  if (overlap_cnt > 0) break;
764  }
765  }
766 // verify_patch_integrity((*p_it).first, &((*p_it).second), patches, info);
767  }
768  }
769  return total_overlapping;
770 }
771 
772 
773 static void
774 bot_partition(struct Manifold_Info *info)
775 {
776  std::map< size_t, std::set<size_t> > face_groups;
777  std::map< size_t, std::set<size_t> > patches;
778  // Calculate face normals dot product with bounding rpp planes
779  for (size_t i=0; i < info->bot->num_faces; ++i) {
780  size_t pt_A, pt_B, pt_C;
781  size_t result_max;
782  double vdot = 0.0;
783  double result = 0.0;
784  vect_t a, b, norm_dir;
785  // Add vert -> edge and edge -> face mappings to global map.
786  pt_A = info->bot->faces[i*3+0]*3;
787  pt_B = info->bot->faces[i*3+1]*3;
788  pt_C = info->bot->faces[i*3+2]*3;
789  info->edge_to_face[mk_edge(pt_A, pt_B)].insert(i);
790  info->edge_to_face[mk_edge(pt_B, pt_C)].insert(i);
791  info->edge_to_face[mk_edge(pt_C, pt_A)].insert(i);
792  // Categorize face
793  VSUB2(a, &info->bot->vertices[info->bot->faces[i*3+1]*3], &info->bot->vertices[info->bot->faces[i*3]*3]);
794  VSUB2(b, &info->bot->vertices[info->bot->faces[i*3+2]*3], &info->bot->vertices[info->bot->faces[i*3]*3]);
795  VCROSS(norm_dir, a, b);
796  VUNITIZE(norm_dir);
797  info->face_normals.Append(ON_3dVector(norm_dir[0], norm_dir[1], norm_dir[2]));
798  (*info).face_normals.At((int)i)->Unitize();
799  for (size_t j=0; j < (size_t)(*info).vectors.Count(); j++) {
800  vdot = ON_DotProduct(*(*info).vectors.At((int)j), *(*info).face_normals.At((int)i));
801  info->norm_results[std::make_pair(i, j)] = vdot;
802  if (vdot > result) {
803  result_max = j;
804  result = vdot;
805  }
806  }
807  face_groups[result_max].insert(i);
808  info->face_to_plane[i]=result_max;
809  }
810  // Sort the groups by their face count - we want to start with the group containing
811  // either the least or the most faces - try most.
812  std::map< size_t, std::set<size_t> >::iterator fg_it;
813  std::vector<size_t> ordered_vects((*info).vectors.Count());
814  std::set<int> face_group_set;
815  for (int i = 0; i < (*info).vectors.Count(); i++) {
816  face_group_set.insert(i);
817  }
818  size_t array_pos = 0;
819  while (!face_group_set.empty()) {
820  size_t curr_max = 0;
821  size_t curr_vect = 0;
822  for (fg_it = face_groups.begin(); fg_it != face_groups.end(); fg_it++) {
823  if ((*fg_it).second.size() > curr_max) {
824  if (face_group_set.find((*fg_it).first) != face_group_set.end()) {
825  curr_max = (*fg_it).second.size();
826  curr_vect = (*fg_it).first;
827  }
828  }
829  }
830  ordered_vects.at(array_pos) = curr_vect;
831  face_group_set.erase(curr_vect);
832  array_pos++;
833  }
834 
835  // All faces must belong to some patch - continue until all faces are processed
836  for (int i = 0; i < 6; i++) {
837  fg_it = face_groups.find(ordered_vects.at(i));
838  std::set<size_t> *faces = &((*fg_it).second);
839  while (!faces->empty()) {
840  info->patch_cnt++;
841  // Find largest remaining face in group
842  double face_size_criteria = 0.0;
843  size_t smallest_face = 0;
844  for (std::set<size_t>::iterator fgl_it = faces->begin(); fgl_it != faces->end(); fgl_it++) {
845  double fa = face_area(info->bot, *fgl_it);
846  if (fa > face_size_criteria) {
847  smallest_face = (*fgl_it);
848  face_size_criteria = fa;
849  }
850  }
851  std::queue<size_t> face_queue;
852  FaceList::iterator f_it;
853  face_queue.push(smallest_face);
854  face_groups[info->face_to_plane[face_queue.front()]].erase(face_queue.front());
855  size_t current_plane = info->face_to_plane[face_queue.front()];
856  //size_t start_face_num = face_queue.front();
857  while (!face_queue.empty()) {
858  size_t face_num = face_queue.front();
859  face_queue.pop();
860  patches[info->patch_cnt].insert(face_num);
861  info->patch_to_plane[info->patch_cnt] = current_plane;
862  info->face_to_patch[face_num] = info->patch_cnt;
863 
864  std::set<size_t> connected_faces;
865  std::set<size_t>::iterator cf_it;
866  get_connected_faces(info->bot, face_num, &(info->edge_to_face), &connected_faces);
867  for (cf_it = connected_faces.begin(); cf_it != connected_faces.end() ; cf_it++) {
868  double curr_face_area = face_area(info->bot, (*cf_it));
869  if (curr_face_area < (face_size_criteria*10) && curr_face_area > (face_size_criteria*0.1)) {
870  if (face_groups[info->face_to_plane[(*cf_it)]].find((*cf_it)) != face_groups[info->face_to_plane[(*cf_it)]].end()) {
871  if (info->face_to_plane[(*cf_it)] == current_plane) {
872  // if (ON_DotProduct( *(*info).face_normals.At((start_face_num)),*(*info).face_normals.At((*cf_it)) ) > 0.5) {
873  // Large patches pose a problem for feature preservation - make an attempt to ensure "large"
874  // patches are flat.
875  if (patches[info->patch_cnt].size() > info->patch_size_threshold) {
876  vect_t origin;
877  size_t ok = 1;
878  VMOVE(origin, &info->bot->vertices[info->bot->faces[(face_num)*3]*3]);
879  ON_Plane plane(ON_3dPoint(origin[0], origin[1], origin[2]), *(*info).face_normals.At((int)(face_num)));
880  for (int pt = 0; pt < 3; pt++) {
881  point_t cpt;
882  VMOVE(cpt, &info->bot->vertices[info->bot->faces[((*cf_it))*3+pt]*3]);
883  double dist_to_plane = plane.DistanceTo(ON_3dPoint(cpt[0], cpt[1], cpt[2]));
884  //std::cout << "Distance[" << face_num << "," << (*cf_it) << "](" << pt << "): " << dist_to_plane << "\n";
885  if (dist_to_plane > 0) {
886  double dist = DIST_PT_PT(origin, cpt);
887  double angle = atan(dist_to_plane/dist);
888  if (angle > info->neighbor_angle_threshold) ok = 0;
889  }
890  }
891  if (ok) {
892  face_queue.push((*cf_it));
893  face_groups[info->face_to_plane[(*cf_it)]].erase((*cf_it));
894  }
895  } else {
896  face_queue.push((*cf_it));
897  face_groups[info->face_to_plane[(*cf_it)]].erase((*cf_it));
898  }
899  // }
900  }
901  }
902  }
903  }
904  }
905  }
906  }
907 
908  // Find "interior" patches and see if they can be merged with the parent patch
909 
910 
911  // "Shave" sharp triangles off of the edges by shifting them from one patch to
912  // another - doesn't avoid all sharp corners, but does produce some cleanup.
913 
914  std::map< size_t, std::set<size_t> >::iterator p_it;
915  size_t shaved_cnt = 1;
916  size_t edge_smoothing_count = 0;
917  while (shaved_cnt != 0 && edge_smoothing_count <= 20) {
918  shaved_cnt = 0;
919  edge_smoothing_count++;
920  for (p_it = patches.begin(); p_it != patches.end(); p_it++) {
921  shaved_cnt += shift_edge_triangles(&patches, (*p_it).first, info);
922  }
923  std::cout << "Triangle shifts in smoothing pass " << edge_smoothing_count << ": " << shaved_cnt << "\n";
924  }
925 
926  // Identify triangles that overlap in the patches associated high-level planar projection, and
927  // assign them to their own patches - they pose a problem for surface fitting if left in the
928  // original overlapping patches.
929  size_t overlapping_face_cnt = 0;
930  overlapping_face_cnt = overlapping_edge_triangles(&patches, info);
931  std::cout << "Found " << overlapping_face_cnt << " overlapping faces in projection\n";
932 
933  // now, populate &(info->patches) with the non-empty patches and fix *_to_patch (maybe make that
934  // part of sync_structure_maps.
935  std::ofstream patch_map;
936  patch_map.open("patch_map.txt");
937  for (p_it = patches.begin(); p_it != patches.end(); p_it++) {
938  if (patches[(*p_it).first].size() > 0) {
939  patch_map << (*p_it).first << " -> " << info->patches.size() << "\n";
940  info->patches[info->patches.size()].insert(patches[(*p_it).first].begin(), patches[(*p_it).first].end());
941  }
942  }
943  patch_map.close();
944 
945  std::cout << "Patch count: " << info->patches.size() << "\n";
946  plot_faces(&(info->patches), info, "patches.pl");
947  plot_patch_borders(&(info->patches), info, "patch_borders.pl");
948 
949  // Now that the mesh is fully partitioned, build the maps needed for curve discovery
950  sync_structure_maps(info);
951 }
952 
953 
954 /**********************************************************************************
955  * Code for moving from patches with edge segments to a structured topology
956  * of curve segments and loops.
957  **********************************************************************************/
958 static void
959 find_edges(struct Manifold_Info *info)
960 {
961  std::map< size_t, size_t> vert_to_m_V; // Keep mapping between vertex index being used in Brep and BoT index
962  std::map< size_t, std::set<size_t> >::iterator p_it;
963  std::map< size_t, EdgeList > all_patch_edges;
964  for (p_it = info->patches.begin(); p_it != info->patches.end(); p_it++) {
965  if (!info->patches[(*p_it).first].empty()) {
966  find_edge_segments(&((*p_it).second), &(all_patch_edges[(*p_it).first]), info);
967  }
968  }
969  FILE* pcurveplot = fopen("polycurves.pl", "w");
970  FILE *nurbs_curves = fopen("nurbs_curves.pl", "w");
971  FILE *curves_plot = fopen("curves.pl", "w");
972  for (p_it = info->patches.begin(); p_it != info->patches.end(); p_it++) {
973  if (!info->patches[(*p_it).first].empty()) {
974  EdgeList *patch_edges = &(all_patch_edges[(*p_it).first]);
975  while (!patch_edges->empty()) {
976  // We know the current patch - find out what the other one is for
977  // this edge.
978  const Edge *first_edge = &(*(patch_edges->begin()));
979  std::set<size_t> edge_patches = (*(info->edge_to_patch.find(*first_edge))).second;
980  edge_patches.erase((*p_it).first);
981  size_t other_patch_id = *(edge_patches.begin());
982  // Now we know our patches - we need to start with the first line
983  // segment, and build up a set of edges until one of the halting
984  // criteria is met.
985  std::set<Edge> polycurve_edges;
986  std::set<Edge>::iterator pc_it;
987  std::queue<Edge> edge_queue;
988  edge_queue.push(*(patch_edges->begin()));
989  all_patch_edges[(*p_it).first].erase(edge_queue.front());
990  all_patch_edges[other_patch_id].erase(edge_queue.front());
991  while (!edge_queue.empty()) {
992  // Add the first edge in the queue
993  Edge curr_edge = edge_queue.front();
994  edge_queue.pop();
995  polycurve_edges.insert(curr_edge);
996  //std::cout << "Current edge: (" << curr_edge.first << "," << curr_edge.second << ")\n";
997 
998  // Using the current edge, identify candidate edges and evaluate them.
999  std::set<Edge> eset;
1000  std::set<Edge>::iterator e_it;
1001  // If a vertex is used by three patches, it is a stopping point for curve
1002  // buildup. Otherwise, consider its edges.
1003  if (info->vert_to_patch[curr_edge.first].size() <= 2) {
1004  eset.insert(info->vert_to_edge[curr_edge.first].begin(), info->vert_to_edge[curr_edge.first].end());
1005  }
1006  if (info->vert_to_patch[curr_edge.second].size() <= 2) {
1007  eset.insert(info->vert_to_edge[curr_edge.second].begin(), info->vert_to_edge[curr_edge.second].end());
1008  }
1009  // We don't need to re-consider any edge we've already considered.
1010  for (e_it = eset.begin(); e_it != eset.end(); e_it++) {
1011  if (polycurve_edges.find((*e_it)) != polycurve_edges.end()) {
1012  eset.erase((*e_it));
1013  }
1014  }
1015  // For the new candidate edges, pull their patches and see if
1016  // they match the current patches. If they do, and the edge has not
1017  // already been removed from one of the patches, this edge is part of
1018  // the current polycurve.
1019  for (e_it = eset.begin(); e_it != eset.end(); e_it++) {
1020  std::set<size_t> *ep = &(info->edge_to_patch[(*e_it)]);
1021  if (ep->find(other_patch_id) != ep->end() && ep->find((*p_it).first) != ep->end() && polycurve_edges.find(*e_it) == polycurve_edges.end()) {
1022  edge_queue.push(*e_it);
1023  all_patch_edges[(*p_it).first].erase(*e_it);
1024  all_patch_edges[other_patch_id].erase(*e_it);
1025  }
1026  }
1027  }
1028  // populate polycurve
1029  std::map<size_t, std::set<size_t> > vert_to_verts;
1030  std::map<size_t, std::set<size_t> >::iterator vv_it;
1031  for (pc_it = polycurve_edges.begin(); pc_it != polycurve_edges.end(); pc_it++) {
1032  vert_to_verts[(*pc_it).first].insert((*pc_it).second);
1033  vert_to_verts[(*pc_it).second].insert((*pc_it).first);
1034  }
1035  // Find start and end points, if any
1036  std::pair<size_t, size_t> start_end;
1037  start_end.first = INT_MAX;
1038  start_end.second = INT_MAX;
1039  for (vv_it = vert_to_verts.begin(); vv_it != vert_to_verts.end(); vv_it++) {
1040  if ((*vv_it).second.size() == 1) {
1041  if (start_end.first == INT_MAX) {
1042  start_end.first = (*vv_it).first;
1043  } else {
1044  if (start_end.second == INT_MAX) {
1045  start_end.second = (*vv_it).first;
1046  }
1047  }
1048  }
1049  }
1050  // Get curve ID number
1051  size_t curve_id = info->polycurves.size();
1052 
1053  // If we have a loop, need different halting conditions for
1054  // curve assembly.
1055  size_t threshold = 0;
1056  if (start_end.first == INT_MAX && start_end.second == INT_MAX) {
1057  threshold = 1;
1058  start_end.first = *(*vert_to_verts.begin()).second.begin();
1059  start_end.second = *(*vert_to_verts.begin()).second.begin();
1060  }
1061 
1062  // Using the starting point, build a vector with the points
1063  // in curve order.
1064  size_t next_pt = start_end.first;
1065  size_t start_end_cnt = 0;
1066  while (start_end_cnt <= threshold) {
1067  if (next_pt == start_end.second) start_end_cnt++;
1068  size_t old_pt = next_pt;
1069  info->polycurves[curve_id].push_back(old_pt);
1070  next_pt = *(vert_to_verts[old_pt].begin());
1071  vert_to_verts[next_pt].erase(old_pt);
1072  }
1073 
1074  // Make m_V vertices out of start and end, m_C3 curve of fit, and BrepEdge using
1075  // start/end pts and curve. When new edge is added, brep->m_E.Count() is the index
1076  // of the new edge. - get it with ON_BrepEdge* edge = brep.Edge(ei); verts are m_vi[0] or m_vi[1],
1077  // or (possibly) edge.Vertex(0)
1078  int vs_id, ve_id;
1079  if (vert_to_m_V.find(start_end.first) == vert_to_m_V.end()) {
1080  ON_BrepVertex& newV = info->brep->NewVertex(ON_3dPoint(&info->bot->vertices[start_end.first]), 0.0001);
1081  vs_id = newV.m_vertex_index;
1082  vert_to_m_V[start_end.first] = vs_id;
1083  } else {
1084  vs_id = vert_to_m_V[start_end.first];
1085  }
1086  if (vert_to_m_V.find(start_end.second) == vert_to_m_V.end()) {
1087  ON_BrepVertex& newV = info->brep->NewVertex(ON_3dPoint(&info->bot->vertices[start_end.second]), 0.0001);
1088  ve_id = newV.m_vertex_index;
1089  vert_to_m_V[start_end.second] = ve_id;
1090  } else {
1091  ve_id = vert_to_m_V[start_end.second];
1092  }
1093  ON_3dPointArray curve_pnts;
1094  std::vector<size_t>::iterator v_it;
1095  // TODO - for curves from patches that had to have triangles removed, need to add
1096  // additional points - a loose fit runs the risk of overlapping NURBS curves in
1097  // projections. May even need a different fit for those situations - worst case,
1098  // use linear edges instead of smooth approximations for all curves that initially failed the
1099  // overlapping triangles test.
1100  for (v_it = info->polycurves[curve_id].begin(); v_it != info->polycurves[curve_id].end(); v_it++) {
1101  ON_3dPoint pt_3d(&info->bot->vertices[(*v_it)]);
1102  curve_pnts.Append(pt_3d);
1103  }
1104  ON_NurbsCurve *curve_nurb;
1105  int c3i;
1106  if (curve_pnts.Count() == 2) {
1107  c3i = info->brep->AddEdgeCurve(new ON_LineCurve(*(curve_pnts.First()), *(curve_pnts.Last())));
1108  } else {
1109  curve_nurb = interpolateLocalCubicCurve(curve_pnts);
1110  c3i = info->brep->AddEdgeCurve(curve_nurb);
1111  }
1112  ON_BrepVertex& StartV = info->brep->m_V[vs_id];
1113  StartV.m_tolerance = 1e-3;
1114  ON_BrepVertex& EndV= info->brep->m_V[ve_id];
1115  EndV.m_tolerance = 1e-3;
1116  ON_BrepEdge& edge = info->brep->NewEdge(StartV, EndV, c3i, NULL ,0);
1117  edge.m_tolerance = 1e-3;
1118 
1119  // Let the patches know they have a curve associated with them.
1120  info->patch_edges[(*p_it).first].insert(edge.m_edge_index);
1121  info->patch_edges[other_patch_id].insert(edge.m_edge_index);
1122 
1123  // Plot curve
1124  int r = int(256*drand48() + 1.0);
1125  int g = int(256*drand48() + 1.0);
1126  int b = int(256*drand48() + 1.0);
1127  plot_curve(info->bot, &(info->polycurves[curve_id]), r, g, b, pcurveplot);
1128  plot_curve(info->bot, &(info->polycurves[curve_id]), r, g, b, curves_plot);
1129 
1130  pl_color(nurbs_curves, r, g, b);
1131  double pt1[3], pt2[3];
1132  int plotres = 50;
1133  const ON_Curve* edge_curve = edge.EdgeCurveOf();
1134  ON_Interval dom = edge_curve->Domain();
1135  for (int i = 1; i <= 50; i++) {
1136  ON_3dPoint p = edge_curve->PointAt(dom.ParameterAt((double)(i - 1)
1137  / (double)plotres));
1138  VMOVE(pt1, p);
1139  p = edge_curve->PointAt(dom.ParameterAt((double) i / (double)plotres));
1140  VMOVE(pt2, p);
1141  pdv_3move(nurbs_curves, pt1);
1142  pdv_3cont(nurbs_curves, pt2);
1143  pdv_3move(curves_plot, pt1);
1144  pdv_3cont(curves_plot, pt2);
1145  }
1146 
1147  }
1148  }
1149  }
1150  fclose(curves_plot);
1151  fclose(nurbs_curves);
1152  fclose(pcurveplot);
1153 }
1154 
1155 
1156 static void
1157 find_outer_loop(std::map<size_t, std::vector<size_t> > *loops, size_t *outer_loop, std::set<size_t> *inner_loops, struct Manifold_Info *info)
1158 {
1159  // If we have more than one loop, we need to identify the outer loop. OpenNURBS handles outer
1160  // and inner loops separately, so we need to know which loop is our outer surface boundary.
1161  if ((*loops).size() == 1) {
1162  *outer_loop = 0;
1163  } else {
1164  std::map<size_t, std::vector<size_t> >::iterator l_it;
1165  fastf_t diag_max = 0.0;
1166  for (l_it = loops->begin(); l_it != loops->end(); l_it++) {
1167  double boxmin[3] = {0.0, 0.0, 0.0};
1168  double boxmax[3] = {0.0, 0.0, 0.0};
1169  std::vector<size_t>::iterator v_it;
1170  for (v_it = (*l_it).second.begin(); v_it != (*l_it).second.end(); v_it++) {
1171  ON_BrepEdge& edge = info->brep->m_E[(int)(*v_it)];
1172  const ON_Curve* edge_curve = edge.EdgeCurveOf();
1173  edge_curve->GetBBox((double *)&boxmin, (double *)&boxmax, 1);
1174  }
1175  ON_3dPoint pmin(boxmin[0], boxmin[1], boxmin[2]);
1176  ON_3dPoint pmax(boxmax[0], boxmax[1], boxmax[2]);
1177  if (pmin.DistanceTo(pmax) > diag_max) {
1178  *outer_loop = (*l_it).first;
1179  diag_max = pmin.DistanceTo(pmax);
1180  }
1181  }
1182  for (l_it = loops->begin(); l_it != loops->end(); l_it++) {
1183  if ((*l_it).first != *outer_loop) {
1184  (*inner_loops).insert((*l_it).first);
1185  }
1186  }
1187  }
1188 }
1189 
1190 
1191 // Do pullbacks to generate trimming curves in 2D, use them to create BrepTrim and BrepLoop structures for
1192 // outer and inner trimming loops
1193 static void
1194 build_loop(size_t patch_id, size_t loop_index, ON_BrepLoop::TYPE loop_type, std::map<size_t, std::vector<size_t> > *loops, struct Manifold_Info *info)
1195 {
1196  ON_BrepFace& face = info->brep->m_F[(int)patch_id];
1197  const ON_Surface *surface = face.SurfaceOf();
1198  ON_Interval xdom = surface->Domain(0);
1199  ON_Interval ydom = surface->Domain(1);
1200  // Start with outer loop
1201  ON_BrepLoop& loop = info->brep->NewLoop(loop_type, face);
1202  // build surface tree
1203  brlcad::SurfaceTree* st = new brlcad::SurfaceTree(&face, false);
1204 
1205  std::vector<size_t>::iterator loop_it;
1206  std::vector<size_t> *loop_edges = &((*loops)[loop_index]);
1207  int vert_prev = -1;
1208  //std::cout << "Patch " << patch_id << " loop edges: \n";
1209  bool trim_rev = false;
1210  for (loop_it = loop_edges->begin(); loop_it != loop_edges->end(); loop_it++) {
1211  size_t curr_edge = (*loop_it);
1212  // Will we need to flip the trim?
1213  ON_BrepEdge& edge = info->brep->m_E[(int)curr_edge];
1214  if (vert_prev != -1) {
1215  if (vert_prev == edge.m_vi[0]) {
1216  trim_rev = false;
1217  } else {
1218  trim_rev = true;
1219  }
1220  }
1221  if (trim_rev) {
1222  vert_prev = edge.m_vi[0];
1223  } else {
1224  vert_prev = edge.m_vi[1];
1225  }
1226  //std::cout << "edge verts: " << edge.m_vi[0] << "," << edge.m_vi[1] << " flip=" << trim_rev << "\n";
1227 
1228  ON_2dPointArray curve_pnts_2d;
1229  const ON_Curve* edge_curve = edge.EdgeCurveOf();
1230  ON_Interval dom = edge_curve->Domain();
1231 
1232  ON_3dPoint pt_3d = edge_curve->PointAt(dom.ParameterAt(0));
1233  ON_2dPoint pt_2d;
1234  ON_2dPoint loop_anchor;
1235  ON_2dPoint pt_2d_prev;
1236  int prev_trim_rev = 0;
1237  int istart = 1;
1238 
1239  size_t pullback_failures;
1240  if (loop_it == loop_edges->begin()) {
1241  pullback_failures = 0;
1242  prev_trim_rev = trim_rev;
1243  int found_first_pt = 0;
1244  while (!found_first_pt && istart < 50) {
1245  pt_3d = edge_curve->PointAt(dom.ParameterAt((double)(istart-1)/(double)50));
1246  ON_3dPoint p3d_pullback = ON_3dPoint::UnsetPoint;
1247  double distance = DBL_MAX;
1248  if (surface_GetClosestPoint3dFirstOrder(face.SurfaceOf(),pt_3d,pt_2d,p3d_pullback,distance,0,BREP_SAME_POINT_TOLERANCE,10/* within tolerance */)) {
1249  if (xdom.Includes(pt_2d.x) && ydom.Includes(pt_2d.y)) {
1250  curve_pnts_2d.Append(pt_2d);
1251  loop_anchor = pt_2d;
1252  found_first_pt = 1;
1253  } else {
1254  istart++;
1255  }
1256  } else {
1257  std::cout << "Pullback failure on first pt (" << patch_id << "," << curr_edge << "," << (double)(istart-1)/(double)50 << "):\n";
1258  std::cout << "Expected: " << pt_3d.x << "," << pt_3d.y << "," << pt_3d.z << "\n";
1259  std::cout << "Got: " << p3d_pullback.x << "," << p3d_pullback.y << "," << p3d_pullback.z << "\n";
1260  std::cout << "Got(2D): " << pt_2d.x << "," << pt_2d.y << "\n";
1261  ON_2dPoint pt_2d_cp;
1262  (void)get_closest_point(pt_2d_cp, &face, pt_3d, st);
1263  std::cout << "get_closest_point(2D): " << pt_2d_cp.x << "," << pt_2d_cp.y << "\n";
1264  pullback_failures++;
1265  istart++;
1266  }
1267  }
1268  } else {
1269  curve_pnts_2d.Append(pt_2d_prev);
1270  }
1271  // XXX todo: dynamically sample the curve - must use consistent method for all sampling, else
1272  // surface may not contain points sought by curve
1273  if (!trim_rev) {
1274  for (int i = istart; i < 50; i++) {
1275  pt_3d = edge_curve->PointAt(dom.ParameterAt((double)(i)/(double)50));
1276  ON_3dPoint p3d_pullback = ON_3dPoint::UnsetPoint;
1277  double distance = DBL_MAX;
1278  if (surface_GetClosestPoint3dFirstOrder(face.SurfaceOf(),pt_3d,pt_2d,p3d_pullback,distance,0,BREP_SAME_POINT_TOLERANCE,10/* within tolerance */)) {
1279  if (xdom.Includes(pt_2d.x) && ydom.Includes(pt_2d.y)) {
1280  curve_pnts_2d.Append(pt_2d);
1281  pt_2d_prev = pt_2d;
1282  }
1283  } else {
1284  std::cout << "Pullback failure (" << patch_id << "," << curr_edge << "," << (double)(i)/(double)50 << "):\n";
1285  std::cout << "Expected: " << pt_3d.x << "," << pt_3d.y << "," << pt_3d.z << "\n";
1286  std::cout << "Got: " << p3d_pullback.x << "," << p3d_pullback.y << "," << p3d_pullback.z << "\n";
1287  std::cout << "Got(2D): " << pt_2d.x << "," << pt_2d.y << "\n";
1288  ON_2dPoint pt_2d_cp;
1289  (void)get_closest_point(pt_2d_cp, &face, pt_3d, st);
1290  std::cout << "get_closest_point(2D): " << pt_2d_cp.x << "," << pt_2d_cp.y << "\n";
1291  pullback_failures++;
1292  }
1293  }
1294  } else {
1295  for (int i = 50; i > istart; i--) {
1296  pt_3d = edge_curve->PointAt(dom.ParameterAt((double)(i)/(double)50));
1297  ON_3dPoint p3d_pullback = ON_3dPoint::UnsetPoint;
1298  double distance = DBL_MAX;
1299  if (surface_GetClosestPoint3dFirstOrder(face.SurfaceOf(),pt_3d,pt_2d,p3d_pullback,distance,0,BREP_SAME_POINT_TOLERANCE,10/* within tolerance */)) {
1300  if (xdom.Includes(pt_2d.x) && ydom.Includes(pt_2d.y)) {
1301  curve_pnts_2d.Append(pt_2d);
1302  pt_2d_prev = pt_2d;
1303  }
1304  } else {
1305  std::cout << "Pullback failure (" << patch_id << "," << curr_edge << "," << (double)(i)/(double)50 << "):\n";
1306  std::cout << "Expected: " << pt_3d.x << "," << pt_3d.y << "," << pt_3d.z << "\n";
1307  std::cout << "Got: " << p3d_pullback.x << "," << p3d_pullback.y << "," << p3d_pullback.z << "\n";
1308  std::cout << "Got(2D): " << pt_2d.x << "," << pt_2d.y << "\n";
1309  ON_2dPoint pt_2d_cp;
1310  (void)get_closest_point(pt_2d_cp, &face, pt_3d, st);
1311  std::cout << "get_closest_point(2D): " << pt_2d_cp.x << "," << pt_2d_cp.y << "\n";
1312  pullback_failures++;
1313  }
1314  }
1315  }
1316  // For final curve, doesn't matter what last pullback is - we MUST force the loop to close.
1317  if (loop_it+1 == loop_edges->end()) {
1318  curve_pnts_2d.Append(loop_anchor);
1319  }
1320  //if (trim_rev) {curve_pnts_2d.Reverse();}
1321  ON_Curve *trim_curve = interpolateCurve(curve_pnts_2d);
1322  int c2i = info->brep->AddTrimCurve(trim_curve);
1323  ON_BrepTrim& trim = info->brep->NewTrim(edge, trim_rev, loop, c2i);
1324  trim.m_type = ON_BrepTrim::mated;
1325  trim.m_tolerance[0] = 1e-3;
1326  trim.m_tolerance[1] = 1e-3;
1327  }
1328  if (info->brep->LoopDirection(loop) != 1 && loop_type == ON_BrepLoop::outer) {
1329  info->brep->FlipLoop(loop);
1330  }
1331  if (info->brep->LoopDirection(loop) != -1 && loop_type == ON_BrepLoop::inner) {
1332  info->brep->FlipLoop(loop);
1333  }
1334 }
1335 
1336 
1337 static void
1338 find_loops(struct Manifold_Info *info)
1339 {
1340  std::map< size_t, std::set<size_t> >::iterator p_it;
1341  for (p_it = info->patch_edges.begin(); p_it != info->patch_edges.end(); p_it++) {
1342  std::map<size_t, std::vector<size_t> > loops;
1343  std::map<size_t, std::set<size_t> > vert_to_edges;
1344  struct bu_vls name;
1345  bu_vls_init(&name);
1346  bu_vls_printf(&name, "loops_patch_%d.pl", (int)(*p_it).first);
1347  FILE* ploopplot = fopen(bu_vls_addr(&name), "w");
1348  std::set<size_t> curr_edges = (*p_it).second;
1349  for (std::set<size_t>::iterator edge_it = curr_edges.begin(); edge_it != curr_edges.end(); edge_it++) {
1350  ON_BrepEdge& edge = info->brep->m_E[(int)(*edge_it)];
1351  if (edge.m_vi[0] == edge.m_vi[1]) {
1352  size_t curr_loop = loops.size();
1353  loops[curr_loop].push_back((*edge_it));
1354  plot_loop(&loops[curr_loop], info, ploopplot);
1355  curr_edges.erase(*edge_it);
1356  } else {
1357  vert_to_edges[edge.m_vi[0]].insert(*edge_it);
1358  vert_to_edges[edge.m_vi[1]].insert(*edge_it);
1359  }
1360  }
1361 
1362  while (curr_edges.size() > 0) {
1363  size_t curr_loop = loops.size();
1364  std::queue<size_t> edge_queue;
1365  edge_queue.push(*(curr_edges.begin()));
1366  curr_edges.erase(edge_queue.front());
1367  int vert_to_match = info->brep->m_E[(int)edge_queue.front()].m_vi[0];
1368  while (!edge_queue.empty()) {
1369  size_t curr_edge = edge_queue.front();
1370  edge_queue.pop();
1371  loops[curr_loop].push_back(curr_edge);
1372  ON_BrepEdge& edge = info->brep->m_E[(int)curr_edge];
1373  if (vert_to_match == edge.m_vi[0]) {
1374  vert_to_match = edge.m_vi[1];
1375  } else {
1376  vert_to_match = edge.m_vi[0];
1377  }
1378  // use vert_to_edges to assemble other curves
1379  std::set<size_t> candidate_edges;
1380  candidate_edges.insert(vert_to_edges[vert_to_match].begin(), vert_to_edges[vert_to_match].end());
1381  candidate_edges.erase(curr_edge);
1382  for (std::set<size_t>::iterator c_it = candidate_edges.begin(); c_it != candidate_edges.end(); c_it++) {
1383  if (curr_edges.find(*c_it) != curr_edges.end()) {
1384  edge_queue.push(*c_it);
1385  curr_edges.erase(*c_it);
1386  }
1387  }
1388  }
1389 #ifdef BOT_TO_NURBS_DEBUG
1390  std::cout << "Patch " << (*p_it).first << " loop " << curr_loop << ": ";
1391  for (size_t i = 0; i < loops[curr_loop].size(); i++) {
1392  std::cout << loops[curr_loop][i] << ",";
1393  }
1394  std::cout << "\n";
1395 #endif
1396  plot_loop(&loops[curr_loop], info, ploopplot);
1397  }
1398  fclose(ploopplot);
1399  // Now that we have loops, determine which is the outer loop
1400  size_t outer_loop;
1401  std::set<size_t> inner_loops;
1402  find_outer_loop(&loops, &outer_loop, &inner_loops, info);
1403  if (loops.size() > 1) {
1404  std::cout << "Patch " << (*p_it).first << " outer loop: " << outer_loop << "\n";
1405  }
1406 
1407  build_loop((*p_it).first, outer_loop, ON_BrepLoop::outer, &loops, info);
1408  // If there are any inner loops, build them too
1409  for (std::set<size_t>::iterator il_it = inner_loops.begin(); il_it != inner_loops.end(); il_it++) {
1410  build_loop((*p_it).first, (*il_it), ON_BrepLoop::inner, &loops, info);
1411  }
1412  }
1413 }
1414 
1415 
1416 /**********************************************************************************
1417  * Code for fitted 3D NURBS surfaces to patches.
1418  **********************************************************************************/
1419 
1420 static void
1421 PatchToVector3d(struct rt_bot_internal *bot, size_t curr_patch, struct Manifold_Info *info, on_fit::vector_vec3d &data)
1422 {
1423  std::set<size_t> *faces = &(info->patches[curr_patch]);
1424  std::set<size_t>::iterator f_it;
1425  std::set<size_t> verts;
1426  unsigned int i = 0;
1427  //ON_3dPointArray pnts;
1428  for (i = 0; i < bot->num_vertices; i++) {
1429  //printf("v(%d): %f %f %f\n", i, V3ARGS(&bot->vertices[3*i]));
1430  }
1431  for (f_it = faces->begin(); f_it != faces->end(); f_it++) {
1432  verts.insert(bot->faces[(*f_it)*3+0]);
1433  verts.insert(bot->faces[(*f_it)*3+1]);
1434  verts.insert(bot->faces[(*f_it)*3+2]);
1435  }
1436  for (f_it = verts.begin(); f_it != verts.end(); f_it++) {
1437  //printf("vert %d\n", (int)(*f_it));
1438  //printf("vert(%d): %f %f %f\n", (int)(*f_it), V3ARGS(&bot->vertices[(*f_it)*3]));
1439  data.push_back(ON_3dVector(V3ARGS(&bot->vertices[(*f_it)*3])));
1440  //pnts.Append(ON_3dPoint(V3ARGS(&bot->vertices[(*f_it)*3])));
1441  }
1442 
1443  // Points on the edges of surfaces make problems for our current 2D uv pullback routine. Uncomment
1444  // the following to try to force the surface edges away from the actual patch boundaries by adding
1445  // "extension" points to the fit.
1446 #if 0
1447  ON_BoundingBox bbox;
1448  pnts.GetTightBoundingBox(bbox);
1449  fastf_t diagonal = bbox.Diagonal().Length();
1450  ON_Plane best_fit_plane;
1451  fit_plane(curr_patch, faces, info, &best_fit_plane);
1452 
1453  fastf_t greatest_distance_to_plane = 0.0;
1454  for (int j = 0; j < pnts.Count(); ++j) {
1455  fastf_t dist = best_fit_plane.DistanceTo(pnts[j]);
1456  if (dist > greatest_distance_to_plane) greatest_distance_to_plane = dist;
1457  }
1458  ON_3dVector move = best_fit_plane.Normal();
1459  double vdot = ON_DotProduct(move, *info->face_normals.At((int)*(info->patches[curr_patch].begin())));
1460  if (vdot > 0.05) move.Reverse();
1461  ON_3dPoint newcenter = best_fit_plane.Origin() + move * greatest_distance_to_plane;
1462  best_fit_plane.SetOrigin(newcenter);
1463 
1464  struct bu_vls filename;
1465  bu_vls_init(&filename);
1466  bu_vls_sprintf(&filename, "extensions_%d.pl", curr_patch);
1467  FILE* edge_plot = fopen(bu_vls_addr(&filename), "w");
1468  EdgeList patch_edges_extend;
1469  EdgeList::iterator e_it;
1470  point_t pt1, pt2;
1471  VSET(pt1, best_fit_plane.Origin().x, best_fit_plane.Origin().y, best_fit_plane.Origin().z)
1472  ON_Circle bounding_circle(best_fit_plane, diagonal);
1473  for (int t = 0; t < 10; t++) {
1474  data.push_back(ON_3dVector(bounding_circle.PointAt(t*36)));
1475  VSET(pt2, bounding_circle.PointAt(t*36).x, bounding_circle.PointAt(t*36).y, bounding_circle.PointAt(t*36).z);
1476 
1477  pdv_3move(edge_plot, pt1);
1478  pdv_3cont(edge_plot, pt2);
1479  }
1480  fclose(edge_plot);
1481 #endif
1482 
1483  // Edges are important for patch merging - tighten them by adding more edge
1484  // points than just the vertex points. Use points from the 3D NURBS edge curves
1485  // instead of the bot edges to ensure the surface includes the volume needed
1486  // for curve pullback.
1487  std::set<size_t> *patch_edges = &(info->patch_edges[curr_patch]);
1488  std::set<size_t>::iterator pe_it;
1489  for (pe_it = patch_edges->begin(); pe_it != patch_edges->end(); pe_it++) {
1490  ON_BrepEdge& edge = info->brep->m_E[(int)(*pe_it)];
1491  const ON_Curve* edge_curve = edge.EdgeCurveOf();
1492  ON_Interval dom = edge_curve->Domain();
1493  // XXX todo: dynamically sample the curve
1494  ON_3dPoint pt_3d = edge_curve->PointAt(dom.ParameterAt(0));
1495  data.push_back(ON_3dVector(pt_3d));
1496  for (int ipts = 1; ipts < 50; ipts++) {
1497  pt_3d = edge_curve->PointAt(dom.ParameterAt((double)(ipts)/(double)50));
1498  data.push_back(ON_3dVector(pt_3d));
1499  }
1500  pt_3d = edge_curve->PointAt(dom.ParameterAt(1.0));
1501  data.push_back(ON_3dVector(pt_3d));
1502  }
1503 }
1504 
1505 
1506 // Actually fit the NURBS surfaces and build faces
1507 static void
1508 find_surfaces(struct Manifold_Info *info)
1509 {
1510  std::map< size_t, std::set<size_t> >::iterator p_it;
1511  for (p_it = info->patches.begin(); p_it != info->patches.end(); p_it++) {
1512  std::cout << "Found patch " << (*p_it).first << " with " << (*p_it).second.size() << " faces\n";
1513  unsigned order(3);
1515  PatchToVector3d(info->bot, (*p_it).first, info, data.interior);
1516  ON_NurbsSurface nurbs = on_fit::FittingSurface::initNurbsPCABoundingBox(order, &data);
1517  on_fit::FittingSurface fit(&data, nurbs);
1519  params.interior_smoothness = 0.15;
1520  params.interior_weight = 1.0;
1521  params.boundary_smoothness = 0.15;
1522  params.boundary_weight = 0.0;
1523  // NURBS refinement
1524  for (unsigned i = 0; i < 5; i++) {
1525  fit.refine(0);
1526  fit.refine(1);
1527  }
1528 
1529  // fitting iterations
1530  for (unsigned i = 0; i < 2; i++) {
1531  fit.assemble(params);
1532  fit.solve();
1533  }
1534  ON_BrepFace *face2 = info->brep2->NewFace(fit.m_nurbs);
1535  info->brep2->m_S.Append(new ON_NurbsSurface(fit.m_nurbs));
1536  int si = info->brep->AddSurface(new ON_NurbsSurface(fit.m_nurbs));
1537  ON_BrepFace &face = info->brep->NewFace(si);
1538  // Check face normal against plane
1539  ON_Interval udom = fit.m_nurbs.Domain(0);
1540  ON_Interval vdom = fit.m_nurbs.Domain(1);
1541  ON_3dVector surf_norm = fit.m_nurbs.NormalAt(udom.ParameterAt(0.5), vdom.ParameterAt(0.5));
1542  double vdot = ON_DotProduct(surf_norm, *info->face_normals.At((int)*(info->patches[(*p_it).first].begin())));
1543  if (vdot < 0) {
1544  info->brep->FlipFace(face);
1545  info->brep2->FlipFace(*face2);
1546  }
1547  }
1548 }
1549 
1550 
1551 int
1552 main(int argc, char *argv[])
1553 {
1554  struct db_i *dbip;
1555  struct directory *dp;
1556  struct rt_db_internal intern;
1557  struct rt_bot_internal *bot_ip = NULL;
1558  struct rt_wdb *wdbp;
1559  struct bu_vls name;
1560  struct bu_vls bname;
1561 
1562  FaceList::iterator f_it;
1563 
1564  Manifold_Info info;
1565  info.vectors.Append(ON_3dVector(-1,0,0));
1566  info.vectors.Append(ON_3dVector(0,-1,0));
1567  info.vectors.Append(ON_3dVector(0,0,-1));
1568  info.vectors.Append(ON_3dVector(1,0,0));
1569  info.vectors.Append(ON_3dVector(0,1,0));
1570  info.vectors.Append(ON_3dVector(0,0,1));
1571  // Smaller values here will do better at preserving detail but will increase patch count.
1572  info.patch_size_threshold = 100;
1573  // Smaller values mean more detail preservation, but result in higher patch counts.
1574  info.neighbor_angle_threshold = 0.05;
1575  // A "normal" value for the "parallel to projection plane" face check should be around 0.55
1576  info.face_plane_parallel_threshold = 0.55;
1577  // Initialize the patch count variable
1578  info.patch_cnt = -1;
1579  bu_vls_init(&name);
1580 
1581  if (argc != 3) {
1582  bu_exit(1, "Usage: %s file.g object", argv[0]);
1583  }
1584 
1585  dbip = db_open(argv[1], DB_OPEN_READWRITE);
1586  if (dbip == DBI_NULL) {
1587  bu_exit(1, "ERROR: Unable to read from %s\n", argv[1]);
1588  }
1589 
1590  if (db_dirbuild(dbip) < 0)
1591  bu_exit(1, "ERROR: Unable to read from %s\n", argv[1]);
1592 
1593  dp = db_lookup(dbip, argv[2], LOOKUP_QUIET);
1594  if (dp == RT_DIR_NULL) {
1595  bu_exit(1, "ERROR: Unable to look up object %s\n", argv[2]);
1596  }
1597 
1598  RT_DB_INTERNAL_INIT(&intern)
1599  if (rt_db_get_internal(&intern, dp, dbip, NULL, &rt_uniresource) < 0) {
1600  bu_exit(1, "ERROR: Unable to get internal representation of %s\n", argv[2]);
1601  }
1602 
1603  if (intern.idb_minor_type != DB5_MINORTYPE_BRLCAD_BOT) {
1604  bu_exit(1, "ERROR: object %s does not appear to be of type BoT\n", argv[2]);
1605  } else {
1606  bot_ip = (struct rt_bot_internal *)intern.idb_ptr;
1607  }
1608  RT_BOT_CK_MAGIC(bot_ip);
1609 
1610  info.bot = bot_ip;
1611 
1612  // Do the initial partitioning of the mesh into patches
1613  bu_log("Performing initial partitioning of %s\n", dp->d_namep);
1614  bot_partition(&info);
1615 
1616  // Create the Brep data structure to hold curves and topology information
1617  info.brep = ON_Brep::New();
1618  info.brep2 = ON_Brep::New();
1619 
1620 
1621  // Now, using the patch sets, construct brep edges
1622  find_edges(&info);
1623 
1624  find_surfaces(&info);
1625 
1626  // Build the loops
1627  find_loops(&info);
1628 
1629  // Generate BRL-CAD brep object
1630  wdbp = wdb_dbopen(dbip, RT_WDB_TYPE_DB_DISK);
1631  bu_vls_init(&bname);
1632  bu_vls_sprintf(&bname, "%s_brep", argv[2]);
1633  if (mk_brep(wdbp, bu_vls_addr(&bname), info.brep) == 0) {
1634  bu_log("Generated brep object %s\n", bu_vls_addr(&bname));
1635  }
1636 
1637  // For debugging, generate individual face BReps
1638  for (int fc = 0; fc < info.brep->m_F.Count(); fc++) {
1639  ON_Brep *brep_face = info.brep->DuplicateFace(fc, false);
1640  bu_vls_sprintf(&bname, "%s_face_%d", argv[2], fc);
1641  if (mk_brep(wdbp, bu_vls_addr(&bname), brep_face) == 0) {
1642  bu_log("Generated brep object %s\n", bu_vls_addr(&bname));
1643  }
1644  }
1645 
1646 
1647 //#ifdef BOT_TO_NURBS_DEBUG
1648  for (std::map< size_t, std::set<size_t> >::iterator np_it = info.patches.begin(); np_it != info.patches.end(); np_it++) {
1649  struct bu_vls pname;
1650  bu_vls_init(&pname);
1651  bu_vls_printf(&pname, "patch_%d.pl", (int)(*np_it).first);
1652  FILE* plot_file = fopen(bu_vls_addr(&pname), "w");
1653  int r = int(256*drand48() + 1.0);
1654  int g = int(256*drand48() + 1.0);
1655  int b = int(256*drand48() + 1.0);
1656  std::set<size_t> *nfaces = &((*np_it).second);
1657  std::set<size_t>::iterator nf_it;
1658  for (nf_it = nfaces->begin(); nf_it != nfaces->end(); nf_it++) {
1659  plot_face(&bot_ip->vertices[bot_ip->faces[(*nf_it)*3+0]*3], &bot_ip->vertices[bot_ip->faces[(*nf_it)*3+1]*3], &bot_ip->vertices[bot_ip->faces[(*nf_it)*3+2]*3], r, g ,b, plot_file);
1660  }
1661  fclose(plot_file);
1662  }
1663 //#endif
1664 
1665 #if 0
1666  FILE* curve_plot = fopen("curve_plot.pl", "w");
1667  pl_color(curve_plot, 255, 255, 255);
1668  std::set<std::pair<size_t, size_t> > patch_interactions;
1669  std::set<std::pair<size_t, size_t> >::iterator o_it;
1670  for (p_it = info.patches.begin(); p_it != info.patches.end(); p_it++) {
1671  EdgeList::iterator e_it;
1672  for (e_it = info.patch_edges[(*p_it).first].begin(); e_it != info.patch_edges[(*p_it).first].end(); e_it++) {
1673  std::set<size_t> ep = info.edge_to_patch[(*e_it)];
1674  ep.erase((*p_it).first);
1675  if ((*p_it).first < (*ep.begin())) {
1676  patch_interactions.insert(std::make_pair((*p_it).first, (*ep.begin())));
1677  } else {
1678  patch_interactions.insert(std::make_pair((*ep.begin()), (*p_it).first));
1679  }
1680  }
1681  }
1682  std::cout << "Patch interaction count: " << patch_interactions.size() << "\n";
1683  for (o_it = patch_interactions.begin(); o_it != patch_interactions.end(); o_it++) {
1684  std::cout << "Intersecting " << (*o_it).first << " with " << (*o_it).second << "\n";
1685  ON_SimpleArray<ON_NurbsCurve*> intersect3d;
1686  ON_SimpleArray<ON_NurbsCurve*> intersect_uv2d;
1687  ON_SimpleArray<ON_NurbsCurve*> intersect_st2d;
1688  brlcad::surface_surface_intersection((*(info.surface_array.At(info.patch_to_surface[(*o_it).first]))), (*(info.surface_array.At(info.patch_to_surface[(*o_it).second]))), intersect3d, intersect_uv2d, intersect_st2d);
1689  for (int k = 0; k < intersect3d.Count(); k++) {
1690  plot_ncurve(*(intersect3d[k]),curve_plot);
1691  delete intersect3d[k];
1692  }
1693 
1694  }
1695  fclose(curve_plot);
1696 #endif
1697  return 0;
1698 }
1699 
1700 // Local Variables:
1701 // tab-width: 8
1702 // mode: C++
1703 // c-basic-offset: 4
1704 // indent-tabs-mode: t
1705 // c-file-style: "stroustrup"
1706 // End:
1707 // ex: shiftwidth=4 tabstop=8
EdgeToFace edge_to_face
void bu_vls_init(struct bu_vls *vp)
Definition: vls.c:56
char * d_namep
pointer to name string
Definition: raytrace.h:859
char filename[MAXLENGTH]
Definition: human.c:105
Definition: raytrace.h:800
std::map< size_t, std::set< size_t > > patches
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
int rt_db_get_internal(struct rt_db_internal *ip, const struct directory *dp, const struct db_i *dbip, const mat_t mat, struct resource *resp)
Definition: dir.c:76
ON_Curve * interpolateCurve(ON_2dPointArray &samples)
int main(int argc, char *argv[])
std::map< size_t, std::set< Edge > > VertToEdge
EdgeToPatch edge_to_patch
FaceToPatch face_to_patch
std::map< size_t, std::set< size_t > > VertToPatch
void pdv_3move(register FILE *plotfp, const fastf_t *pt)
Definition: plot3.c:618
static ON_NurbsSurface initNurbsPCABoundingBox(int order, NurbsDataSurface *data, ON_3dVector z=ON_3dVector(0.0, 0.0, 1.0))
Initializing a B-Spline surface using principal-component-analysis and bounding box of points...
Definition: clone.c:90
double neighbor_angle_threshold
std::pair< size_t, size_t > Edge
Fitting a B-Spline surface to 3D point-clouds using point-distance-minimization Based on paper: TODO...
size_t patch_size_threshold
#define VSET(a, b, c, d)
Definition: color.c:53
ON_3dVectorArray vectors
struct directory * db_lookup(const struct db_i *, const char *name, int noisy)
Definition: db_lookup.c:153
#define st
Header file for the BRL-CAD common definitions.
std::map< size_t, size_t > FaceToPatch
std::map< std::pair< size_t, size_t >, fastf_t > norm_results
std::set< Edge > EdgeList
bool get_closest_point(ON_2dPoint &outpt, ON_BrepFace *face, const ON_3dPoint &point, SurfaceTree *tree, double tolerance)
std::map< size_t, size_t > face_to_plane
std::pair< size_t, size_t > mk_edge(size_t pt_A, size_t pt_B)
#define DB_OPEN_READWRITE
Definition: raytrace.h:3555
COMPLEX data[64]
Definition: fftest.c:34
vector_vec3d interior
< input
struct resource rt_uniresource
default. Defined in librt/globals.c
Definition: globals.c:41
ON_NurbsCurve * interpolateLocalCubicCurve(ON_2dPointArray &Q)
void bu_exit(int status, const char *fmt,...) _BU_ATTR_NORETURN _BU_ATTR_PRINTF23
Definition: bomb.c:195
ON_3dVectorArray face_normals
double face_plane_parallel_threshold
#define RT_DB_INTERNAL_INIT(_p)
Definition: raytrace.h:199
void bu_vls_sprintf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:707
#define LOOKUP_QUIET
Definition: raytrace.h:893
Parameters for fitting.
#define V3ARGS(a)
Definition: color.c:56
#define RT_WDB_TYPE_DB_DISK
Definition: raytrace.h:1295
VertToEdge vert_to_edge
Data structure for NURBS surface fitting (FittingSurface, FittingSurfaceTDM, FittingCylinder, GlobalOptimization, GlobalOptimizationTDM)
Definition: opennurbs_fit.h:96
std::vector< ON_3dVector > vector_vec3d
Definition: opennurbs_fit.h:92
int mk_brep(struct rt_wdb *wdbp, const char *name, ON_Brep *brep)
Definition: brep.cpp:42
void plot_face(ON_3dPoint *pt1, ON_3dPoint *pt2, ON_3dPoint *pt3, int r, int g, int b, FILE *c_plot)
void pdv_3cont(register FILE *plotfp, const fastf_t *pt)
Definition: plot3.c:630
#define UNUSED(parameter)
Definition: common.h:239
std::set< size_t > FaceList
std::map< size_t, size_t > PatchToPlane
char * bu_vls_addr(const struct bu_vls *vp)
Definition: vls.c:111
double drand48(void)
void pl_color(register FILE *plotfp, int r, int g, int b)
Definition: plot3.c:325
bool surface_GetClosestPoint3dFirstOrder(const ON_Surface *surf, const ON_3dPoint &p, ON_2dPoint &p2d, ON_3dPoint &p3d, double &current_distance, int quadrant, double same_point_tol, double within_distance_tol)
std::map< size_t, std::vector< size_t > > polycurves
std::map< size_t, std::set< size_t > > patch_edges
PatchToPlane patch_to_plane
std::map< Edge, std::set< size_t > > EdgeToPatch
void * idb_ptr
Definition: raytrace.h:195
int overlap(struct application *ap, struct partition *pp, struct region *reg1, struct region *reg2, struct partition *hp)
Definition: gqa.c:843
int bn_tri_tri_isect_coplanar(point_t V0, point_t V1, point_t V2, point_t U0, point_t U1, point_t U2, int area_flag)
Tomas Möller's triangle/triangle intersection routines from the article.
Definition: tri_tri.c:182
#define DBI_NULL
Definition: raytrace.h:827
std::map< Edge, std::set< size_t > > EdgeToFace
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
#define RT_DIR_NULL
Definition: raytrace.h:875
struct db_i * db_open(const char *name, const char *mode)
Definition: db_open.c:59
int idb_minor_type
ID_xxx.
Definition: raytrace.h:193
#define A
Definition: msr.c:51
int db_dirbuild(struct db_i *dbip)
Definition: db5_scan.c:301
struct rt_bot_internal * bot
Definition: vls.h:56
double fastf_t
Definition: defines.h:300
struct rt_wdb * wdb_dbopen(struct db_i *dbip, int mode)
Definition: wdb.c:64
VertToPatch vert_to_patch