Blender  V3.3
octree.cpp
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
3 #ifdef WITH_CXX_GUARDEDALLOC
4 # include "MEM_guardedalloc.h"
5 #endif
6 
7 #include "octree.h"
8 #include <Eigen/Dense>
9 #include <limits>
10 #include <time.h>
11 
18 /* set to non-zero value to enable debugging output */
19 #define DC_DEBUG 0
20 
21 #if DC_DEBUG
22 /* enable debug printfs */
23 # define dc_printf printf
24 #else
25 /* disable debug printfs */
26 # define dc_printf(...) \
27  do { \
28  } while (0)
29 #endif
30 
32  DualConAllocOutput alloc_output_func,
33  DualConAddVert add_vert_func,
34  DualConAddQuad add_quad_func,
35  DualConFlags flags,
36  DualConMode dualcon_mode,
37  int depth,
38  float threshold,
39  float sharpness)
40  : use_flood_fill(flags & DUALCON_FLOOD_FILL),
41  /* note on `use_manifold':
42 
43  After playing around with this option, the only case I could
44  find where this option gives different results is on
45  relatively thin corners. Sometimes along these corners two
46  vertices from separate sides will be placed in the same
47  position, so hole gets filled with a 5-sided face, where two
48  of those vertices are in the same 3D location. If
49  `use_manifold' is disabled, then the modifier doesn't
50  generate two separate vertices so the results end up as all
51  quads.
52 
53  Since the results are just as good with all quads, there
54  doesn't seem any reason to allow this to be toggled in
55  Blender. -nicholasbishop
56  */
57  use_manifold(false),
58  hermite_num(sharpness),
59  mode(dualcon_mode),
60  alloc_output(alloc_output_func),
61  add_vert(add_vert_func),
62  add_quad(add_quad_func)
63 {
64  thresh = threshold;
65  reader = mr;
66  dimen = 1 << GRID_DIMENSION;
68  nodeCount = nodeSpace = 0;
69  maxDepth = depth;
70  mindimen = (dimen >> maxDepth);
72  buildTable();
73 
75 
76  // Initialize memory
77 #ifdef IN_VERBOSE_MODE
78  dc_printf("Range: %f origin: %f, %f,%f \n", range, origin[0], origin[1], origin[2]);
79  dc_printf("Initialize memory...\n");
80 #endif
81  initMemory();
82  root = (Node *)createInternal(0);
83 
84  // Read MC table
85 #ifdef IN_VERBOSE_MODE
86  dc_printf("Reading contour table...\n");
87 #endif
88  cubes = new Cubes();
89 }
90 
92 {
93  delete cubes;
94  freeMemory();
95 }
96 
98 {
99  // Scan triangles
100 #if DC_DEBUG
101  clock_t start, finish;
102  start = clock();
103 #endif
104 
105  addAllTriangles();
106  resetMinimalEdges();
107  preparePrimalEdgesMask(&root->internal);
108 
109 #if DC_DEBUG
110  finish = clock();
111  dc_printf("Time taken: %f seconds \n",
112  (double)(finish - start) / CLOCKS_PER_SEC);
113 #endif
114 
115  // Generate signs
116  // Find holes
117 #if DC_DEBUG
118  dc_printf("Patching...\n");
119  start = clock();
120 #endif
121  trace();
122 #if DC_DEBUG
123  finish = clock();
124  dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
125 # ifdef IN_VERBOSE_MODE
126  dc_printf("Holes: %d Average Length: %f Max Length: %d \n",
127  numRings,
128  (float)totRingLengths / (float)numRings,
129  maxRingLength);
130 # endif
131 #endif
132 
133  // Check again
134  int tnumRings = numRings;
135  trace();
136 #ifdef IN_VERBOSE_MODE
137  dc_printf("Holes after patching: %d \n", numRings);
138 #endif
139  numRings = tnumRings;
140 
141 #if DC_DEBUG
142  dc_printf("Building signs...\n");
143  start = clock();
144 #endif
145  buildSigns();
146 #if DC_DEBUG
147  finish = clock();
148  dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
149 #endif
150 
151  if (use_flood_fill) {
152  /*
153  start = clock();
154  floodFill();
155  // Check again
156  tnumRings = numRings;
157  trace();
158  dc_printf("Holes after filling: %d \n", numRings);
159  numRings = tnumRings;
160  buildSigns();
161  finish = clock();
162  dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
163  */
164 #if DC_DEBUG
165  start = clock();
166  dc_printf("Removing components...\n");
167 #endif
168  floodFill();
169  buildSigns();
170  // dc_printf("Checking...\n");
171  // floodFill();
172 #if DC_DEBUG
173  finish = clock();
174  dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
175 #endif
176  }
177 
178  // Output
179 #if DC_DEBUG
180  start = clock();
181 #endif
182  writeOut();
183 #if DC_DEBUG
184  finish = clock();
185 #endif
186  // dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
187 
188  // Print info
189 #ifdef IN_VERBOSE_MODE
190  printMemUsage();
191 #endif
192 }
193 
194 void Octree::initMemory()
195 {
200 
210 }
211 
212 void Octree::freeMemory()
213 {
214  for (int i = 0; i < 9; i++) {
215  alloc[i]->destroy();
216  delete alloc[i];
217  }
218 
219  for (int i = 0; i < 4; i++) {
220  leafalloc[i]->destroy();
221  delete leafalloc[i];
222  }
223 }
224 
225 void Octree::printMemUsage()
226 {
227  int totalbytes = 0;
228  dc_printf("********* Internal nodes: \n");
229  for (int i = 0; i < 9; i++) {
230  alloc[i]->printInfo();
231 
232  totalbytes += alloc[i]->getAll() * alloc[i]->getBytes();
233  }
234  dc_printf("********* Leaf nodes: \n");
235  int totalLeafs = 0;
236  for (int i = 0; i < 4; i++) {
237  leafalloc[i]->printInfo();
238 
239  totalbytes += leafalloc[i]->getAll() * leafalloc[i]->getBytes();
240  totalLeafs += leafalloc[i]->getAllocated();
241  }
242 
243  dc_printf("Total allocated bytes on disk: %d \n", totalbytes);
244  dc_printf("Total leaf nodes: %d\n", totalLeafs);
245 
246  /* Unused when not debuggining. */
247  (void)totalbytes;
248  (void)totalLeafs;
249 }
250 
251 void Octree::resetMinimalEdges()
252 {
253  cellProcParity(root, 0, maxDepth);
254 }
255 
256 void Octree::addAllTriangles()
257 {
258  Triangle *trian;
259  int count = 0;
260 
261 #if DC_DEBUG
262  int total = reader->getNumTriangles();
263  int unitcount = 1000;
264  dc_printf("\nScan converting to depth %d...\n", maxDepth);
265 #endif
266 
267  srand(0);
268 
269  while ((trian = reader->getNextTriangle()) != NULL) {
270  // Drop triangles
271  {
272  addTriangle(trian, count);
273  }
274  delete trian;
275 
276  count++;
277 
278 #if DC_DEBUG
279  if (count % unitcount == 0) {
280  putchar(13);
281 
282  switch ((count / unitcount) % 4) {
283  case 0:
284  dc_printf("-");
285  break;
286  case 1:
287  dc_printf("/");
288  break;
289  case 2:
290  dc_printf("|");
291  break;
292  case 3:
293  dc_printf("\\");
294  break;
295  }
296 
297  float percent = (float)count / total;
298 
299  /*
300  int totbars = 50;
301  int bars =(int)(percent * totbars);
302  for(int i = 0; i < bars; i ++) {
303  putchar(219);
304  }
305  for(i = bars; i < totbars; i ++) {
306  putchar(176);
307  }
308  */
309 
310  dc_printf(" %d triangles: ", count);
311  dc_printf(" %f%% complete.", 100 * percent);
312  }
313 #endif
314  }
315  putchar(13);
316 }
317 
318 /* Prepare a triangle for insertion into the octree; call the other
319  addTriangle() to (recursively) build the octree */
320 void Octree::addTriangle(Triangle *trian, int triind)
321 {
322  int i, j;
323 
324  /* Project the triangle's coordinates into the grid */
325  for (i = 0; i < 3; i++) {
326  for (j = 0; j < 3; j++)
327  trian->vt[i][j] = dimen * (trian->vt[i][j] - origin[j]) / range;
328  }
329 
330  /* Generate projections */
331  int64_t cube[2][3] = {{0, 0, 0}, {dimen, dimen, dimen}};
332  int64_t trig[3][3];
333  for (i = 0; i < 3; i++) {
334  for (j = 0; j < 3; j++)
335  trig[i][j] = (int64_t)(trian->vt[i][j]);
336  }
337 
338  /* Add triangle to the octree */
339  int64_t errorvec = (int64_t)(0);
340  CubeTriangleIsect *proj = new CubeTriangleIsect(cube, trig, errorvec, triind);
341  root = (Node *)addTriangle(&root->internal, proj, maxDepth);
342 
343  delete proj->inherit;
344  delete proj;
345 }
346 
347 #if 0
348 static void print_depth(int height, int maxDepth)
349 {
350  for (int i = 0; i < maxDepth - height; i++)
351  printf(" ");
352 }
353 #endif
354 
355 InternalNode *Octree::addTriangle(InternalNode *node, CubeTriangleIsect *p, int height)
356 {
357  int i;
358  const int vertdiff[8][3] = {
359  {0, 0, 0}, {0, 0, 1}, {0, 1, -1}, {0, 0, 1}, {1, -1, -1}, {0, 0, 1}, {0, 1, -1}, {0, 0, 1}};
360  unsigned char boxmask = p->getBoxMask();
361  CubeTriangleIsect *subp = new CubeTriangleIsect(p);
362 
363  int count = 0;
364  int tempdiff[3] = {0, 0, 0};
365 
366  /* Check triangle against each of the input node's children */
367  for (i = 0; i < 8; i++) {
368  tempdiff[0] += vertdiff[i][0];
369  tempdiff[1] += vertdiff[i][1];
370  tempdiff[2] += vertdiff[i][2];
371 
372  /* Quick pruning using bounding box */
373  if (boxmask & (1 << i)) {
374  subp->shift(tempdiff);
375  tempdiff[0] = tempdiff[1] = tempdiff[2] = 0;
376 
377  /* Pruning using intersection test */
378  if (subp->isIntersecting()) {
379  if (!node->has_child(i)) {
380  if (height == 1)
381  node = addLeafChild(node, i, count, createLeaf(0));
382  else
383  node = addInternalChild(node, i, count, createInternal(0));
384  }
385  Node *chd = node->get_child(count);
386 
387  if (node->is_child_leaf(i))
388  node->set_child(count, (Node *)updateCell(&chd->leaf, subp));
389  else
390  node->set_child(count, (Node *)addTriangle(&chd->internal, subp, height - 1));
391  }
392  }
393 
394  if (node->has_child(i))
395  count++;
396  }
397 
398  delete subp;
399 
400  return node;
401 }
402 
403 LeafNode *Octree::updateCell(LeafNode *node, CubeTriangleIsect *p)
404 {
405  int i;
406 
407  // Edge connectivity
408  int mask[3] = {0, 4, 8};
409  int oldc = 0, newc = 0;
410  float offs[3];
411  float a[3], b[3], c[3];
412 
413  for (i = 0; i < 3; i++) {
414  if (!getEdgeParity(node, mask[i])) {
415  if (p->isIntersectingPrimary(i)) {
416  // actualQuads ++;
417  setEdge(node, mask[i]);
418  offs[newc] = p->getIntersectionPrimary(i);
419  a[newc] = (float)p->inherit->norm[0];
420  b[newc] = (float)p->inherit->norm[1];
421  c[newc] = (float)p->inherit->norm[2];
422  newc++;
423  }
424  }
425  else {
426  offs[newc] = getEdgeOffsetNormal(node, oldc, a[newc], b[newc], c[newc]);
427 
428  oldc++;
429  newc++;
430  }
431  }
432 
433  if (newc > oldc) {
434  // New offsets added, update this node
435  node = updateEdgeOffsetsNormals(node, oldc, newc, offs, a, b, c);
436  }
437 
438  return node;
439 }
440 
441 void Octree::preparePrimalEdgesMask(InternalNode *node)
442 {
443  int count = 0;
444  for (int i = 0; i < 8; i++) {
445  if (node->has_child(i)) {
446  if (node->is_child_leaf(i))
447  createPrimalEdgesMask(&node->get_child(count)->leaf);
448  else
449  preparePrimalEdgesMask(&node->get_child(count)->internal);
450 
451  count++;
452  }
453  }
454 }
455 
456 void Octree::trace()
457 {
458  int st[3] = {
459  0,
460  0,
461  0,
462  };
463  numRings = 0;
464  totRingLengths = 0;
465  maxRingLength = 0;
466 
467  PathList *chdpath = NULL;
468  root = trace(root, st, dimen, maxDepth, chdpath);
469 
470  if (chdpath != NULL) {
471  dc_printf("there are incomplete rings.\n");
472  printPaths(chdpath);
473  }
474 }
475 
476 Node *Octree::trace(Node *newnode, int *st, int len, int depth, PathList *&paths)
477 {
478  len >>= 1;
479  PathList *chdpaths[8];
480  Node *chd[8];
481  int nst[8][3];
482  int i, j;
483 
484  // Get children paths
485  int chdleaf[8];
486  newnode->internal.fill_children(chd, chdleaf);
487 
488  // int count = 0;
489  for (i = 0; i < 8; i++) {
490  for (j = 0; j < 3; j++) {
491  nst[i][j] = st[j] + len * vertmap[i][j];
492  }
493 
494  if (chd[i] == NULL || newnode->internal.is_child_leaf(i)) {
495  chdpaths[i] = NULL;
496  }
497  else {
498  trace(chd[i], nst[i], len, depth - 1, chdpaths[i]);
499  }
500  }
501 
502  // Get connectors on the faces
503  PathList *conn[12];
504  Node *nf[2];
505  int lf[2];
506  int df[2] = {depth - 1, depth - 1};
507  int *nstf[2];
508 
509  newnode->internal.fill_children(chd, chdleaf);
510 
511  for (i = 0; i < 12; i++) {
512  int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
513 
514  for (int j = 0; j < 2; j++) {
515  lf[j] = chdleaf[c[j]];
516  nf[j] = chd[c[j]];
517  nstf[j] = nst[c[j]];
518  }
519 
520  conn[i] = NULL;
521 
522  findPaths((Node **)nf, lf, df, nstf, depth - 1, cellProcFaceMask[i][2], conn[i]);
523 
524 #if 0
525  if (conn[i]) {
526  printPath(conn[i]);
527  }
528 #endif
529  }
530 
531  // Connect paths
532  PathList *rings = NULL;
533  combinePaths(chdpaths[0], chdpaths[1], conn[8], rings);
534  combinePaths(chdpaths[2], chdpaths[3], conn[9], rings);
535  combinePaths(chdpaths[4], chdpaths[5], conn[10], rings);
536  combinePaths(chdpaths[6], chdpaths[7], conn[11], rings);
537 
538  combinePaths(chdpaths[0], chdpaths[2], conn[4], rings);
539  combinePaths(chdpaths[4], chdpaths[6], conn[5], rings);
540  combinePaths(chdpaths[0], NULL, conn[6], rings);
541  combinePaths(chdpaths[4], NULL, conn[7], rings);
542 
543  combinePaths(chdpaths[0], chdpaths[4], conn[0], rings);
544  combinePaths(chdpaths[0], NULL, conn[1], rings);
545  combinePaths(chdpaths[0], NULL, conn[2], rings);
546  combinePaths(chdpaths[0], NULL, conn[3], rings);
547 
548  // By now, only chdpaths[0] and rings have contents
549 
550  // Process rings
551  if (rings) {
552  // printPath(rings);
553 
554  /* Let's count first */
555  PathList *trings = rings;
556  while (trings) {
557  numRings++;
558  totRingLengths += trings->length;
559  if (trings->length > maxRingLength) {
560  maxRingLength = trings->length;
561  }
562  trings = trings->next;
563  }
564 
565  // printPath(rings);
566  newnode = patch(newnode, st, (len << 1), rings);
567  }
568 
569  // Return incomplete paths
570  paths = chdpaths[0];
571  return newnode;
572 }
573 
574 void Octree::findPaths(
575  Node *node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList *&paths)
576 {
577  if (!(node[0] && node[1])) {
578  return;
579  }
580 
581  if (!(leaf[0] && leaf[1])) {
582  // Not at the bottom, recur
583 
584  // Fill children nodes
585  int i, j;
586  Node *chd[2][8];
587  int chdleaf[2][8];
588  int nst[2][8][3];
589 
590  for (j = 0; j < 2; j++) {
591  if (!leaf[j]) {
592  node[j]->internal.fill_children(chd[j], chdleaf[j]);
593 
594  int len = (dimen >> (maxDepth - depth[j] + 1));
595  for (i = 0; i < 8; i++) {
596  for (int k = 0; k < 3; k++) {
597  nst[j][i][k] = st[j][k] + len * vertmap[i][k];
598  }
599  }
600  }
601  }
602 
603  // 4 face calls
604  Node *nf[2];
605  int df[2];
606  int lf[2];
607  int *nstf[2];
608  for (i = 0; i < 4; i++) {
609  int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
610  for (int j = 0; j < 2; j++) {
611  if (leaf[j]) {
612  lf[j] = leaf[j];
613  nf[j] = node[j];
614  df[j] = depth[j];
615  nstf[j] = st[j];
616  }
617  else {
618  lf[j] = chdleaf[j][c[j]];
619  nf[j] = chd[j][c[j]];
620  df[j] = depth[j] - 1;
621  nstf[j] = nst[j][c[j]];
622  }
623  }
624  findPaths(nf, lf, df, nstf, maxdep - 1, faceProcFaceMask[dir][i][2], paths);
625  }
626  }
627  else {
628  // At the bottom, check this face
629  int ind = (depth[0] == maxdep ? 0 : 1);
630  int fcind = 2 * dir + (1 - ind);
631  if (getFaceParity((LeafNode *)node[ind], fcind)) {
632  // Add into path
633  PathElement *ele1 = new PathElement;
634  PathElement *ele2 = new PathElement;
635 
636  ele1->pos[0] = st[0][0];
637  ele1->pos[1] = st[0][1];
638  ele1->pos[2] = st[0][2];
639 
640  ele2->pos[0] = st[1][0];
641  ele2->pos[1] = st[1][1];
642  ele2->pos[2] = st[1][2];
643 
644  ele1->next = ele2;
645  ele2->next = NULL;
646 
647  PathList *lst = new PathList;
648  lst->head = ele1;
649  lst->tail = ele2;
650  lst->length = 2;
651  lst->next = paths;
652  paths = lst;
653 
654  // int l =(dimen >> maxDepth);
655  }
656  }
657 }
658 
659 void Octree::combinePaths(PathList *&list1, PathList *list2, PathList *paths, PathList *&rings)
660 {
661  // Make new list of paths
662  PathList *nlist = NULL;
663 
664  // Search for each connectors in paths
665  PathList *tpaths = paths;
666  PathList *tlist, *pre;
667  while (tpaths) {
668  PathList *singlist = tpaths;
669  PathList *templist;
670  tpaths = tpaths->next;
671  singlist->next = NULL;
672 
673  // Look for hookup in list1
674  tlist = list1;
675  pre = NULL;
676  while (tlist) {
677  if ((templist = combineSinglePath(list1, pre, tlist, singlist, NULL, singlist)) != NULL) {
678  singlist = templist;
679  continue;
680  }
681  pre = tlist;
682  tlist = tlist->next;
683  }
684 
685  // Look for hookup in list2
686  tlist = list2;
687  pre = NULL;
688  while (tlist) {
689  if ((templist = combineSinglePath(list2, pre, tlist, singlist, NULL, singlist)) != NULL) {
690  singlist = templist;
691  continue;
692  }
693  pre = tlist;
694  tlist = tlist->next;
695  }
696 
697  // Look for hookup in nlist
698  tlist = nlist;
699  pre = NULL;
700  while (tlist) {
701  if ((templist = combineSinglePath(nlist, pre, tlist, singlist, NULL, singlist)) != NULL) {
702  singlist = templist;
703  continue;
704  }
705  pre = tlist;
706  tlist = tlist->next;
707  }
708 
709  // Add to nlist or rings
710  if (isEqual(singlist->head, singlist->tail)) {
711  PathElement *temp = singlist->head;
712  singlist->head = temp->next;
713  delete temp;
714  singlist->length--;
715  singlist->tail->next = singlist->head;
716 
717  singlist->next = rings;
718  rings = singlist;
719  }
720  else {
721  singlist->next = nlist;
722  nlist = singlist;
723  }
724  }
725 
726  // Append list2 and nlist to the end of list1
727  tlist = list1;
728  if (tlist != NULL) {
729  while (tlist->next != NULL) {
730  tlist = tlist->next;
731  }
732  tlist->next = list2;
733  }
734  else {
735  tlist = list2;
736  list1 = list2;
737  }
738 
739  if (tlist != NULL) {
740  while (tlist->next != NULL) {
741  tlist = tlist->next;
742  }
743  tlist->next = nlist;
744  }
745  else {
746  tlist = nlist;
747  list1 = nlist;
748  }
749 }
750 
751 PathList *Octree::combineSinglePath(PathList *&head1,
752  PathList *pre1,
753  PathList *&list1,
754  PathList *&head2,
755  PathList *pre2,
756  PathList *&list2)
757 {
758  if (isEqual(list1->head, list2->head) || isEqual(list1->tail, list2->tail)) {
759  // Reverse the list
760  if (list1->length < list2->length) {
761  // Reverse list1
762  PathElement *prev = list1->head;
764  prev->next = NULL;
765  while (next != NULL) {
766  PathElement *tnext = next->next;
767  next->next = prev;
768 
769  prev = next;
770  next = tnext;
771  }
772 
773  list1->tail = list1->head;
774  list1->head = prev;
775  }
776  else {
777  // Reverse list2
778  PathElement *prev = list2->head;
780  prev->next = NULL;
781  while (next != NULL) {
782  PathElement *tnext = next->next;
783  next->next = prev;
784 
785  prev = next;
786  next = tnext;
787  }
788 
789  list2->tail = list2->head;
790  list2->head = prev;
791  }
792  }
793 
794  if (isEqual(list1->head, list2->tail)) {
795 
796  // Easy case
797  PathElement *temp = list1->head->next;
798  delete list1->head;
799  list2->tail->next = temp;
800 
801  PathList *nlist = new PathList;
802  nlist->length = list1->length + list2->length - 1;
803  nlist->head = list2->head;
804  nlist->tail = list1->tail;
805  nlist->next = NULL;
806 
807  deletePath(head1, pre1, list1);
808  deletePath(head2, pre2, list2);
809 
810  return nlist;
811  }
812  else if (isEqual(list1->tail, list2->head)) {
813  // Easy case
814  PathElement *temp = list2->head->next;
815  delete list2->head;
816  list1->tail->next = temp;
817 
818  PathList *nlist = new PathList;
819  nlist->length = list1->length + list2->length - 1;
820  nlist->head = list1->head;
821  nlist->tail = list2->tail;
822  nlist->next = NULL;
823 
824  deletePath(head1, pre1, list1);
825  deletePath(head2, pre2, list2);
826 
827  return nlist;
828  }
829 
830  return NULL;
831 }
832 
833 void Octree::deletePath(PathList *&head, PathList *pre, PathList *&curr)
834 {
835  PathList *temp = curr;
836  curr = temp->next;
837  delete temp;
838 
839  if (pre == NULL) {
840  head = curr;
841  }
842  else {
843  pre->next = curr;
844  }
845 }
846 
847 void Octree::printElement(PathElement *ele)
848 {
849  if (ele != NULL) {
850  dc_printf("(%d %d %d)", ele->pos[0], ele->pos[1], ele->pos[2]);
851  }
852 }
853 
854 void Octree::printPath(PathList *path)
855 {
856  PathElement *n = path->head;
857  int same = 0;
858 
859 #if DC_DEBUG
860  int len = (dimen >> maxDepth);
861 #endif
862  while (n && (same == 0 || n != path->head)) {
863  same++;
864  dc_printf("(%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len);
865  n = n->next;
866  }
867 
868  if (n == path->head) {
869  dc_printf(" Ring!\n");
870  }
871  else {
872  dc_printf(" %p end!\n", n);
873  }
874 }
875 
876 void Octree::printPath(PathElement *path)
877 {
878  PathElement *n = path;
879  int same = 0;
880 #if DC_DEBUG
881  int len = (dimen >> maxDepth);
882 #endif
883  while (n && (same == 0 || n != path)) {
884  same++;
885  dc_printf("(%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len);
886  n = n->next;
887  }
888 
889  if (n == path) {
890  dc_printf(" Ring!\n");
891  }
892  else {
893  dc_printf(" %p end!\n", n);
894  }
895 }
896 
897 void Octree::printPaths(PathList *path)
898 {
899  PathList *iter = path;
900  int i = 0;
901  while (iter != NULL) {
902  dc_printf("Path %d:\n", i);
903  printPath(iter);
904  iter = iter->next;
905  i++;
906  }
907 }
908 
909 Node *Octree::patch(Node *newnode, int st[3], int len, PathList *rings)
910 {
911 #ifdef IN_DEBUG_MODE
912  dc_printf("Call to PATCH with rings: \n");
913  printPaths(rings);
914 #endif
915 
916  /* Do nothing but couting
917  PathList* tlist = rings;
918  PathList* ttlist;
919  PathElement* telem, * ttelem;
920  while(tlist!= NULL) {
921  // printPath(tlist);
922  numRings ++;
923  totRingLengths += tlist->length;
924  if(tlist->length > maxRingLength) {
925  maxRingLength = tlist->length;
926  }
927  ttlist = tlist;
928  tlist = tlist->next;
929  }
930  return node;
931  */
932 
933  /* Pass onto separate calls in each direction */
934  if (len == mindimen) {
935  dc_printf("Error! should have no list by now.\n");
936  exit(0);
937  }
938 
939  // YZ plane
940  PathList *xlists[2];
941  newnode = patchSplit(newnode, st, len, rings, 0, xlists[0], xlists[1]);
942 
943  // XZ plane
944  PathList *ylists[4];
945  newnode = patchSplit(newnode, st, len, xlists[0], 1, ylists[0], ylists[1]);
946  newnode = patchSplit(newnode, st, len, xlists[1], 1, ylists[2], ylists[3]);
947 
948  // XY plane
949  PathList *zlists[8];
950  newnode = patchSplit(newnode, st, len, ylists[0], 2, zlists[0], zlists[1]);
951  newnode = patchSplit(newnode, st, len, ylists[1], 2, zlists[2], zlists[3]);
952  newnode = patchSplit(newnode, st, len, ylists[2], 2, zlists[4], zlists[5]);
953  newnode = patchSplit(newnode, st, len, ylists[3], 2, zlists[6], zlists[7]);
954 
955  // Recur
956  len >>= 1;
957  int count = 0;
958  for (int i = 0; i < 8; i++) {
959  if (zlists[i] != NULL) {
960  int nori[3] = {
961  st[0] + len * vertmap[i][0], st[1] + len * vertmap[i][1], st[2] + len * vertmap[i][2]};
962  patch(newnode->internal.get_child(count), nori, len, zlists[i]);
963  }
964 
965  if (newnode->internal.has_child(i)) {
966  count++;
967  }
968  }
969 #ifdef IN_DEBUG_MODE
970  dc_printf("Return from PATCH\n");
971 #endif
972  return newnode;
973 }
974 
975 Node *Octree::patchSplit(Node *newnode,
976  int st[3],
977  int len,
978  PathList *rings,
979  int dir,
980  PathList *&nrings1,
981  PathList *&nrings2)
982 {
983 #ifdef IN_DEBUG_MODE
984  dc_printf("Call to PATCHSPLIT with direction %d and rings: \n", dir);
985  printPaths(rings);
986 #endif
987 
988  nrings1 = NULL;
989  nrings2 = NULL;
990  PathList *tmp;
991  while (rings != NULL) {
992  // Process this ring
993  newnode = patchSplitSingle(newnode, st, len, rings->head, dir, nrings1, nrings2);
994 
995  // Delete this ring from the group
996  tmp = rings;
997  rings = rings->next;
998  delete tmp;
999  }
1000 
1001 #ifdef IN_DEBUG_MODE
1002  dc_printf("Return from PATCHSPLIT with \n");
1003  dc_printf("Rings gourp 1:\n");
1004  printPaths(nrings1);
1005  dc_printf("Rings group 2:\n");
1006  printPaths(nrings2);
1007 #endif
1008 
1009  return newnode;
1010 }
1011 
1012 Node *Octree::patchSplitSingle(Node *newnode,
1013  int st[3],
1014  int len,
1015  PathElement *head,
1016  int dir,
1017  PathList *&nrings1,
1018  PathList *&nrings2)
1019 {
1020 #ifdef IN_DEBUG_MODE
1021  dc_printf("Call to PATCHSPLITSINGLE with direction %d and path: \n", dir);
1022  printPath(head);
1023 #endif
1024 
1025  if (head == NULL) {
1026 #ifdef IN_DEBUG_MODE
1027  dc_printf("Return from PATCHSPLITSINGLE with head==NULL.\n");
1028 #endif
1029  return newnode;
1030  }
1031  else {
1032  // printPath(head);
1033  }
1034 
1035  // Walk along the ring to find pair of intersections
1036  PathElement *pre1 = NULL;
1037  PathElement *pre2 = NULL;
1038  int side = findPair(head, st[dir] + len / 2, dir, pre1, pre2);
1039 
1040 #if 0
1041  if (pre1 == pre2) {
1042  int edgelen = (dimen >> maxDepth);
1043  dc_printf("Location: %d %d %d Direction: %d Reso: %d\n",
1044  st[0] / edgelen,
1045  st[1] / edgelen,
1046  st[2] / edgelen,
1047  dir,
1048  len / edgelen);
1049  printPath(head);
1050  exit(0);
1051  }
1052 #endif
1053 
1054  if (side) {
1055  // Entirely on one side
1056  PathList *nring = new PathList();
1057  nring->head = head;
1058 
1059  if (side == -1) {
1060  nring->next = nrings1;
1061  nrings1 = nring;
1062  }
1063  else {
1064  nring->next = nrings2;
1065  nrings2 = nring;
1066  }
1067  }
1068  else {
1069  // Break into two parts
1070  PathElement *nxt1 = pre1->next;
1071  PathElement *nxt2 = pre2->next;
1072  pre1->next = nxt2;
1073  pre2->next = nxt1;
1074 
1075  newnode = connectFace(newnode, st, len, dir, pre1, pre2);
1076 
1077  if (isEqual(pre1, pre1->next)) {
1078  if (pre1 == pre1->next) {
1079  delete pre1;
1080  pre1 = NULL;
1081  }
1082  else {
1083  PathElement *temp = pre1->next;
1084  pre1->next = temp->next;
1085  delete temp;
1086  }
1087  }
1088  if (isEqual(pre2, pre2->next)) {
1089  if (pre2 == pre2->next) {
1090  delete pre2;
1091  pre2 = NULL;
1092  }
1093  else {
1094  PathElement *temp = pre2->next;
1095  pre2->next = temp->next;
1096  delete temp;
1097  }
1098  }
1099 
1100  compressRing(pre1);
1101  compressRing(pre2);
1102 
1103  // Recur
1104  newnode = patchSplitSingle(newnode, st, len, pre1, dir, nrings1, nrings2);
1105  newnode = patchSplitSingle(newnode, st, len, pre2, dir, nrings1, nrings2);
1106  }
1107 
1108 #ifdef IN_DEBUG_MODE
1109  dc_printf("Return from PATCHSPLITSINGLE with \n");
1110  dc_printf("Rings gourp 1:\n");
1111  printPaths(nrings1);
1112  dc_printf("Rings group 2:\n");
1113  printPaths(nrings2);
1114 #endif
1115 
1116  return newnode;
1117 }
1118 
1119 Node *Octree::connectFace(
1120  Node *newnode, int st[3], int len, int dir, PathElement *f1, PathElement *f2)
1121 {
1122 #ifdef IN_DEBUG_MODE
1123  dc_printf("Call to CONNECTFACE with direction %d and length %d path: \n", dir, len);
1124  dc_printf("Path(low side): \n");
1125  printPath(f1);
1126  // checkPath(f1);
1127  dc_printf("Path(high side): \n");
1128  printPath(f2);
1129 // checkPath(f2);
1130 #endif
1131 
1132  // Setup 2D
1133  int pos = st[dir] + len / 2;
1134  int xdir = (dir + 1) % 3;
1135  int ydir = (dir + 2) % 3;
1136 
1137  // Use existing intersections on f1 and f2
1138  int x1, y1, x2, y2;
1139  float p1, q1, p2, q2;
1140 
1141  getFacePoint(f2->next, dir, x1, y1, p1, q1);
1142  getFacePoint(f2, dir, x2, y2, p2, q2);
1143 
1144  float dx = x2 + p2 - x1 - p1;
1145  float dy = y2 + q2 - y1 - q1;
1146 
1147  // Do adapted Bresenham line drawing
1148  float rx = p1, ry = q1;
1149  int incx = 1, incy = 1;
1150  int lx = x1, ly = y1;
1151  int hx = x2, hy = y2;
1152  int choice;
1153  if (x2 < x1) {
1154  incx = -1;
1155  rx = 1 - rx;
1156  lx = x2;
1157  hx = x1;
1158  }
1159  if (y2 < y1) {
1160  incy = -1;
1161  ry = 1 - ry;
1162  ly = y2;
1163  hy = y1;
1164  }
1165 
1166  float sx = dx * incx;
1167  float sy = dy * incy;
1168 
1169  int ori[3];
1170  ori[dir] = pos / mindimen;
1171  ori[xdir] = x1;
1172  ori[ydir] = y1;
1173  int walkdir;
1174  int inc;
1175  float alpha;
1176 
1177  PathElement *curEleN = f1;
1178  PathElement *curEleP = f2->next;
1179  Node *nodeN = NULL, *nodeP = NULL;
1180  LeafNode *curN = locateLeaf(&newnode->internal, len, f1->pos);
1181  LeafNode *curP = locateLeaf(&newnode->internal, len, f2->next->pos);
1182  if (curN == NULL || curP == NULL) {
1183  exit(0);
1184  }
1185  int stN[3], stP[3];
1186  int lenN, lenP;
1187 
1188  /* Unused code, leaving for posterity
1189 
1190  float stpt[3], edpt[3];
1191  stpt[dir] = edpt[dir] =(float) pos;
1192  stpt[xdir] =(x1 + p1) * mindimen;
1193  stpt[ydir] =(y1 + q1) * mindimen;
1194  edpt[xdir] =(x2 + p2) * mindimen;
1195  edpt[ydir] =(y2 + q2) * mindimen;
1196  */
1197  while (ori[xdir] != x2 || ori[ydir] != y2) {
1198  int next;
1199  if (sy * (1 - rx) > sx * (1 - ry)) {
1200  choice = 1;
1201  next = ori[ydir] + incy;
1202  if (next < ly || next > hy) {
1203  choice = 4;
1204  next = ori[xdir] + incx;
1205  }
1206  }
1207  else {
1208  choice = 2;
1209  next = ori[xdir] + incx;
1210  if (next < lx || next > hx) {
1211  choice = 3;
1212  next = ori[ydir] + incy;
1213  }
1214  }
1215 
1216  if (choice & 1) {
1217  ori[ydir] = next;
1218  if (choice == 1) {
1219  rx += (sy == 0 ? 0 : (1 - ry) * sx / sy);
1220  ry = 0;
1221  }
1222 
1223  walkdir = 2;
1224  inc = incy;
1225  alpha = x2 < x1 ? 1 - rx : rx;
1226  }
1227  else {
1228  ori[xdir] = next;
1229  if (choice == 2) {
1230  ry += (sx == 0 ? 0 : (1 - rx) * sy / sx);
1231  rx = 0;
1232  }
1233 
1234  walkdir = 1;
1235  inc = incx;
1236  alpha = y2 < y1 ? 1 - ry : ry;
1237  }
1238 
1239  // Get the exact location of the marcher
1240  int nori[3] = {ori[0] * mindimen, ori[1] * mindimen, ori[2] * mindimen};
1241  float spt[3] = {(float)nori[0], (float)nori[1], (float)nori[2]};
1242  spt[(dir + (3 - walkdir)) % 3] += alpha * mindimen;
1243  if (inc < 0) {
1244  spt[(dir + walkdir) % 3] += mindimen;
1245  }
1246 
1247 #if 0
1248  dc_printf("new x,y: %d %d\n", ori[xdir] / edgelen, ori[ydir] / edgelen);
1249  dc_printf("nori: %d %d %d alpha: %f walkdir: %d\n", nori[0], nori[1], nori[2], alpha, walkdir);
1250  dc_printf("%f %f %f\n", spt[0], spt[1], spt[2]);
1251 #endif
1252 
1253  // Locate the current cells on both sides
1254  newnode = locateCell(&newnode->internal, st, len, nori, dir, 1, nodeN, stN, lenN);
1255  newnode = locateCell(&newnode->internal, st, len, nori, dir, 0, nodeP, stP, lenP);
1256 
1257  updateParent(&newnode->internal, len, st);
1258 
1259  int flag = 0;
1260  // Add the cells to the rings and fill in the patch
1261  PathElement *newEleN;
1262  if (curEleN->pos[0] != stN[0] || curEleN->pos[1] != stN[1] || curEleN->pos[2] != stN[2]) {
1263  if (curEleN->next->pos[0] != stN[0] || curEleN->next->pos[1] != stN[1] ||
1264  curEleN->next->pos[2] != stN[2]) {
1265  newEleN = new PathElement;
1266  newEleN->next = curEleN->next;
1267  newEleN->pos[0] = stN[0];
1268  newEleN->pos[1] = stN[1];
1269  newEleN->pos[2] = stN[2];
1270 
1271  curEleN->next = newEleN;
1272  }
1273  else {
1274  newEleN = curEleN->next;
1275  }
1276  curN = patchAdjacent(&newnode->internal,
1277  len,
1278  curEleN->pos,
1279  curN,
1280  newEleN->pos,
1281  (LeafNode *)nodeN,
1282  walkdir,
1283  inc,
1284  dir,
1285  1,
1286  alpha);
1287 
1288  curEleN = newEleN;
1289  flag++;
1290  }
1291 
1292  PathElement *newEleP;
1293  if (curEleP->pos[0] != stP[0] || curEleP->pos[1] != stP[1] || curEleP->pos[2] != stP[2]) {
1294  if (f2->pos[0] != stP[0] || f2->pos[1] != stP[1] || f2->pos[2] != stP[2]) {
1295  newEleP = new PathElement;
1296  newEleP->next = curEleP;
1297  newEleP->pos[0] = stP[0];
1298  newEleP->pos[1] = stP[1];
1299  newEleP->pos[2] = stP[2];
1300 
1301  f2->next = newEleP;
1302  }
1303  else {
1304  newEleP = f2;
1305  }
1306  curP = patchAdjacent(&newnode->internal,
1307  len,
1308  curEleP->pos,
1309  curP,
1310  newEleP->pos,
1311  (LeafNode *)nodeP,
1312  walkdir,
1313  inc,
1314  dir,
1315  0,
1316  alpha);
1317 
1318  curEleP = newEleP;
1319  flag++;
1320  }
1321 
1322  /*
1323  if(flag == 0) {
1324  dc_printf("error: non-synchronized patching! at \n");
1325  }
1326  */
1327  }
1328 
1329 #ifdef IN_DEBUG_MODE
1330  dc_printf("Return from CONNECTFACE with \n");
1331  dc_printf("Path(low side):\n");
1332  printPath(f1);
1333  checkPath(f1);
1334  dc_printf("Path(high side):\n");
1335  printPath(f2);
1336  checkPath(f2);
1337 #endif
1338 
1339  return newnode;
1340 }
1341 
1342 LeafNode *Octree::patchAdjacent(InternalNode *node,
1343  int len,
1344  int st1[3],
1345  LeafNode *leaf1,
1346  int st2[3],
1347  LeafNode *leaf2,
1348  int walkdir,
1349  int inc,
1350  int dir,
1351  int side,
1352  float alpha)
1353 {
1354 #ifdef IN_DEBUG_MODE
1355  dc_printf("Before patching.\n");
1356  printInfo(st1);
1357  printInfo(st2);
1358  dc_printf(
1359  "-----------------%d %d %d; %d %d %d\n", st1[0], st2[1], st1[2], st2[0], st2[1], st2[2]);
1360 #endif
1361 
1362  // Get edge index on each leaf
1363  int edgedir = (dir + (3 - walkdir)) % 3;
1364  int incdir = (dir + walkdir) % 3;
1365  int ind1 = (edgedir == 1 ? (dir + 3 - edgedir) % 3 - 1 : 2 - (dir + 3 - edgedir) % 3);
1366  int ind2 = (edgedir == 1 ? (incdir + 3 - edgedir) % 3 - 1 : 2 - (incdir + 3 - edgedir) % 3);
1367 
1368  int eind1 = ((edgedir << 2) | (side << ind1) | ((inc > 0 ? 1 : 0) << ind2));
1369  int eind2 = ((edgedir << 2) | (side << ind1) | ((inc > 0 ? 0 : 1) << ind2));
1370 
1371 #ifdef IN_DEBUG_MODE
1372  dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha);
1373  /*
1374  if(alpha < 0 || alpha > 1) {
1375  dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha);
1376  printInfo(st1);
1377  printInfo(st2);
1378  }
1379  */
1380 #endif
1381 
1382  // Flip edge parity
1383  LeafNode *nleaf1 = flipEdge(leaf1, eind1, alpha);
1384  LeafNode *nleaf2 = flipEdge(leaf2, eind2, alpha);
1385 
1386  // Update parent link
1387  updateParent(node, len, st1, nleaf1);
1388  updateParent(node, len, st2, nleaf2);
1389  // updateParent(nleaf1, mindimen, st1);
1390  // updateParent(nleaf2, mindimen, st2);
1391 
1392  /*
1393  float m[3];
1394  dc_printf("Adding new point: %f %f %f\n", spt[0], spt[1], spt[2]);
1395  getMinimizer(leaf1, m);
1396  dc_printf("Cell %d now has minimizer %f %f %f\n", leaf1, m[0], m[1], m[2]);
1397  getMinimizer(leaf2, m);
1398  dc_printf("Cell %d now has minimizer %f %f %f\n", leaf2, m[0], m[1], m[2]);
1399  */
1400 
1401 #ifdef IN_DEBUG_MODE
1402  dc_printf("After patching.\n");
1403  printInfo(st1);
1404  printInfo(st2);
1405 #endif
1406  return nleaf2;
1407 }
1408 
1409 Node *Octree::locateCell(InternalNode *node,
1410  int st[3],
1411  int len,
1412  int ori[3],
1413  int dir,
1414  int side,
1415  Node *&rleaf,
1416  int rst[3],
1417  int &rlen)
1418 {
1419 #ifdef IN_DEBUG_MODE
1420  // dc_printf("Call to LOCATECELL with node ");
1421  // printNode(node);
1422 #endif
1423 
1424  int i;
1425  len >>= 1;
1426  int ind = 0;
1427  for (i = 0; i < 3; i++) {
1428  ind <<= 1;
1429  if (i == dir && side == 1) {
1430  ind |= (ori[i] <= (st[i] + len) ? 0 : 1);
1431  }
1432  else {
1433  ind |= (ori[i] < (st[i] + len) ? 0 : 1);
1434  }
1435  }
1436 
1437 #ifdef IN_DEBUG_MODE
1438 # if 0
1439  dc_printf(
1440  "In LOCATECELL index of ori(%d %d %d) with dir %d side %d in st(%d %d %d, %d) is: %d\n",
1441  ori[0],
1442  ori[1],
1443  ori[2],
1444  dir,
1445  side,
1446  st[0],
1447  st[1],
1448  st[2],
1449  len,
1450  ind);
1451 # endif
1452 #endif
1453 
1454  rst[0] = st[0] + vertmap[ind][0] * len;
1455  rst[1] = st[1] + vertmap[ind][1] * len;
1456  rst[2] = st[2] + vertmap[ind][2] * len;
1457 
1458  if (node->has_child(ind)) {
1459  int count = node->get_child_count(ind);
1460  Node *chd = node->get_child(count);
1461  if (node->is_child_leaf(ind)) {
1462  rleaf = chd;
1463  rlen = len;
1464  }
1465  else {
1466  // Recur
1467  node->set_child(count,
1468  locateCell(&chd->internal, rst, len, ori, dir, side, rleaf, rst, rlen));
1469  }
1470  }
1471  else {
1472  // Create a new child here
1473  if (len == mindimen) {
1474  LeafNode *chd = createLeaf(0);
1475  node = addChild(node, ind, (Node *)chd, 1);
1476  rleaf = (Node *)chd;
1477  rlen = len;
1478  }
1479  else {
1480  // Subdivide the empty cube
1481  InternalNode *chd = createInternal(0);
1482  node = addChild(node, ind, locateCell(chd, rst, len, ori, dir, side, rleaf, rst, rlen), 0);
1483  }
1484  }
1485 
1486 #ifdef IN_DEBUG_MODE
1487  // dc_printf("Return from LOCATECELL with node ");
1488  // printNode(newnode);
1489 #endif
1490  return (Node *)node;
1491 }
1492 
1493 void Octree::checkElement(PathElement * /*ele*/)
1494 {
1495 #if 0
1496  if (ele != NULL && locateLeafCheck(ele->pos) != ele->node) {
1497  dc_printf("Screwed! at pos: %d %d %d\n",
1498  ele->pos[0] >> minshift,
1499  ele->pos[1] >> minshift,
1500  ele->pos[2] >> minshift);
1501  exit(0);
1502  }
1503 #endif
1504 }
1505 
1506 void Octree::checkPath(PathElement *path)
1507 {
1508  PathElement *n = path;
1509  int same = 0;
1510  while (n && (same == 0 || n != path)) {
1511  same++;
1512  checkElement(n);
1513  n = n->next;
1514  }
1515 }
1516 
1517 void Octree::testFacePoint(PathElement *e1, PathElement *e2)
1518 {
1519  int i;
1520  PathElement *e = NULL;
1521  for (i = 0; i < 3; i++) {
1522  if (e1->pos[i] != e2->pos[i]) {
1523  if (e1->pos[i] < e2->pos[i]) {
1524  e = e2;
1525  }
1526  else {
1527  e = e1;
1528  }
1529  break;
1530  }
1531  }
1532 
1533  int x, y;
1534  float p, q;
1535  dc_printf("Test.");
1536  getFacePoint(e, i, x, y, p, q);
1537 }
1538 
1539 void Octree::getFacePoint(PathElement *leaf, int dir, int &x, int &y, float &p, float &q)
1540 {
1541  // Find average intersections
1542  float avg[3] = {0, 0, 0};
1543  float off[3];
1544  int num = 0, num2 = 0;
1545 
1546  LeafNode *leafnode = locateLeaf(leaf->pos);
1547  for (int i = 0; i < 4; i++) {
1548  int edgeind = faceMap[dir * 2][i];
1549  int nst[3];
1550  for (int j = 0; j < 3; j++) {
1551  nst[j] = leaf->pos[j] + mindimen * vertmap[edgemap[edgeind][0]][j];
1552  }
1553 
1554  if (getEdgeIntersectionByIndex(nst, edgeind / 4, off, 1)) {
1555  avg[0] += off[0];
1556  avg[1] += off[1];
1557  avg[2] += off[2];
1558  num++;
1559  }
1560  if (getEdgeParity(leafnode, edgeind)) {
1561  num2++;
1562  }
1563  }
1564  if (num == 0) {
1565  dc_printf("Wrong! dir: %d pos: %d %d %d num: %d\n",
1566  dir,
1567  leaf->pos[0] >> minshift,
1568  leaf->pos[1] >> minshift,
1569  leaf->pos[2] >> minshift,
1570  num2);
1571  avg[0] = (float)leaf->pos[0];
1572  avg[1] = (float)leaf->pos[1];
1573  avg[2] = (float)leaf->pos[2];
1574  }
1575  else {
1576 
1577  avg[0] /= num;
1578  avg[1] /= num;
1579  avg[2] /= num;
1580 
1581  // avg[0] =(float) leaf->pos[0];
1582  // avg[1] =(float) leaf->pos[1];
1583  // avg[2] =(float) leaf->pos[2];
1584  }
1585 
1586  int xdir = (dir + 1) % 3;
1587  int ydir = (dir + 2) % 3;
1588 
1589  float xf = avg[xdir];
1590  float yf = avg[ydir];
1591 
1592 #ifdef IN_DEBUG_MODE
1593  // Is it outside?
1594  // PathElement* leaf = leaf1->len < leaf2->len ? leaf1 : leaf2;
1595 # if 0
1596  float *m = (leaf == leaf1 ? m1 : m2);
1597  if (xf < leaf->pos[xdir] || yf < leaf->pos[ydir] || xf > leaf->pos[xdir] + leaf->len ||
1598  yf > leaf->pos[ydir] + leaf->len) {
1599  dc_printf("Outside cube(%d %d %d), %d : %d %d %f %f\n",
1600  leaf->pos[0],
1601  leaf->pos[1],
1602  leaf->pos[2],
1603  leaf->len,
1604  pos,
1605  dir,
1606  xf,
1607  yf);
1608 
1609  // For now, snap to cell
1610  xf = m[xdir];
1611  yf = m[ydir];
1612  }
1613 # endif
1614 
1615 # if 0
1616  if (alpha < 0 || alpha > 1 || xf < leaf->pos[xdir] || xf > leaf->pos[xdir] + leaf->len ||
1617  yf < leaf->pos[ydir] || yf > leaf->pos[ydir] + leaf->len) {
1618  dc_printf("Alpha: %f Address: %d and %d\n", alpha, leaf1->node, leaf2->node);
1619  dc_printf("GETFACEPOINT result:(%d %d %d) %d min:(%f %f %f);(%d %d %d) %d min:(%f %f %f).\n",
1620  leaf1->pos[0],
1621  leaf1->pos[1],
1622  leaf1->pos[2],
1623  leaf1->len,
1624  m1[0],
1625  m1[1],
1626  m1[2],
1627  leaf2->pos[0],
1628  leaf2->pos[1],
1629  leaf2->pos[2],
1630  leaf2->len,
1631  m2[0],
1632  m2[1],
1633  m2[2]);
1634  dc_printf("Face point at dir %d pos %d: %f %f\n", dir, pos, xf, yf);
1635  }
1636 # endif
1637 #endif
1638 
1639  // Get the integer and float part
1640  x = ((leaf->pos[xdir]) >> minshift);
1641  y = ((leaf->pos[ydir]) >> minshift);
1642 
1643  p = (xf - leaf->pos[xdir]) / mindimen;
1644  q = (yf - leaf->pos[ydir]) / mindimen;
1645 
1646 #ifdef IN_DEBUG_MODE
1647  dc_printf("Face point at dir %d : %f %f\n", dir, xf, yf);
1648 #endif
1649 }
1650 
1651 int Octree::findPair(PathElement *head, int pos, int dir, PathElement *&pre1, PathElement *&pre2)
1652 {
1653  int side = getSide(head, pos, dir);
1654  PathElement *cur = head;
1655  PathElement *anchor;
1656  PathElement *ppre1, *ppre2;
1657 
1658  // Start from this face, find a pair
1659  anchor = cur;
1660  ppre1 = cur;
1661  cur = cur->next;
1662  while (cur != anchor && (getSide(cur, pos, dir) == side)) {
1663  ppre1 = cur;
1664  cur = cur->next;
1665  }
1666  if (cur == anchor) {
1667  // No pair found
1668  return side;
1669  }
1670 
1671  side = getSide(cur, pos, dir);
1672  ppre2 = cur;
1673  cur = cur->next;
1674  while (getSide(cur, pos, dir) == side) {
1675  ppre2 = cur;
1676  cur = cur->next;
1677  }
1678 
1679  // Switch pre1 and pre2 if we start from the higher side
1680  if (side == -1) {
1681  cur = ppre1;
1682  ppre1 = ppre2;
1683  ppre2 = cur;
1684  }
1685 
1686  pre1 = ppre1;
1687  pre2 = ppre2;
1688 
1689  return 0;
1690 }
1691 
1692 int Octree::getSide(PathElement *e, int pos, int dir)
1693 {
1694  return (e->pos[dir] < pos ? -1 : 1);
1695 }
1696 
1697 int Octree::isEqual(PathElement *e1, PathElement *e2)
1698 {
1699  return (e1->pos[0] == e2->pos[0] && e1->pos[1] == e2->pos[1] && e1->pos[2] == e2->pos[2]);
1700 }
1701 
1702 void Octree::compressRing(PathElement *&ring)
1703 {
1704  if (ring == NULL) {
1705  return;
1706  }
1707 #ifdef IN_DEBUG_MODE
1708  dc_printf("Call to COMPRESSRING with path: \n");
1709  printPath(ring);
1710 #endif
1711 
1712  PathElement *cur = ring->next->next;
1713  PathElement *pre = ring->next;
1714  PathElement *prepre = ring;
1715  PathElement *anchor = prepre;
1716 
1717  do {
1718  while (isEqual(cur, prepre)) {
1719  // Delete
1720  if (cur == prepre) {
1721  // The ring has shrinked to a point
1722  delete pre;
1723  delete cur;
1724  anchor = NULL;
1725  break;
1726  }
1727  else {
1728  prepre->next = cur->next;
1729  delete pre;
1730  delete cur;
1731  pre = prepre->next;
1732  cur = pre->next;
1733  anchor = prepre;
1734  }
1735  }
1736 
1737  if (anchor == NULL) {
1738  break;
1739  }
1740 
1741  prepre = pre;
1742  pre = cur;
1743  cur = cur->next;
1744  } while (prepre != anchor);
1745 
1746  ring = anchor;
1747 
1748 #ifdef IN_DEBUG_MODE
1749  dc_printf("Return from COMPRESSRING with path: \n");
1750  printPath(ring);
1751 #endif
1752 }
1753 
1754 void Octree::buildSigns()
1755 {
1756  // First build a lookup table
1757  // dc_printf("Building up look up table...\n");
1758  int size = 1 << 12;
1759  unsigned char table[1 << 12];
1760  for (int i = 0; i < size; i++) {
1761  table[i] = 0;
1762  }
1763  for (int i = 0; i < 256; i++) {
1764  int ind = 0;
1765  for (int j = 11; j >= 0; j--) {
1766  ind <<= 1;
1767  if (((i >> edgemap[j][0]) & 1) ^ ((i >> edgemap[j][1]) & 1)) {
1768  ind |= 1;
1769  }
1770  }
1771 
1772  table[ind] = i;
1773  }
1774 
1775  // Next, traverse the grid
1776  int sg = 1;
1777  int cube[8];
1778  buildSigns(table, root, 0, sg, cube);
1779 }
1780 
1781 void Octree::buildSigns(unsigned char table[], Node *node, int isLeaf, int sg, int rvalue[8])
1782 {
1783  if (node == NULL) {
1784  for (int i = 0; i < 8; i++) {
1785  rvalue[i] = sg;
1786  }
1787  return;
1788  }
1789 
1790  if (isLeaf == 0) {
1791  // Internal node
1792  Node *chd[8];
1793  int leaf[8];
1794  node->internal.fill_children(chd, leaf);
1795 
1796  // Get the signs at the corners of the first cube
1797  rvalue[0] = sg;
1798  int oris[8];
1799  buildSigns(table, chd[0], leaf[0], sg, oris);
1800 
1801  // Get the rest
1802  int cube[8];
1803  for (int i = 1; i < 8; i++) {
1804  buildSigns(table, chd[i], leaf[i], oris[i], cube);
1805  rvalue[i] = cube[i];
1806  }
1807  }
1808  else {
1809  // Leaf node
1810  generateSigns(&node->leaf, table, sg);
1811 
1812  for (int i = 0; i < 8; i++) {
1813  rvalue[i] = getSign(&node->leaf, i);
1814  }
1815  }
1816 }
1817 
1818 void Octree::floodFill()
1819 {
1820  // int threshold =(int)((dimen/mindimen) *(dimen/mindimen) * 0.5f);
1821  int st[3] = {0, 0, 0};
1822 
1823  // First, check for largest component
1824  // size stored in -threshold
1825  clearProcessBits(root, maxDepth);
1826  int threshold = floodFill(root, st, dimen, maxDepth, 0);
1827 
1828  // Next remove
1829  dc_printf("Largest component: %d\n", threshold);
1830  threshold *= thresh;
1831  dc_printf("Removing all components smaller than %d\n", threshold);
1832 
1833  int st2[3] = {0, 0, 0};
1834  clearProcessBits(root, maxDepth);
1835  floodFill(root, st2, dimen, maxDepth, threshold);
1836 }
1837 
1838 void Octree::clearProcessBits(Node *node, int height)
1839 {
1840  int i;
1841 
1842  if (height == 0) {
1843  // Leaf cell,
1844  for (i = 0; i < 12; i++) {
1845  setOutProcess(&node->leaf, i);
1846  }
1847  }
1848  else {
1849  // Internal cell, recur
1850  int count = 0;
1851  for (i = 0; i < 8; i++) {
1852  if (node->internal.has_child(i)) {
1853  clearProcessBits(node->internal.get_child(count), height - 1);
1854  count++;
1855  }
1856  }
1857  }
1858 }
1859 
1860 int Octree::floodFill(LeafNode *leaf, int st[3], int len, int /*height*/, int threshold)
1861 {
1862  int i, j;
1863  int maxtotal = 0;
1864 
1865  // Leaf cell,
1866  int par, inp;
1867 
1868  // Test if the leaf has intersection edges
1869  for (i = 0; i < 12; i++) {
1870  par = getEdgeParity(leaf, i);
1871  inp = isInProcess(leaf, i);
1872 
1873  if (par == 1 && inp == 0) {
1874  // Intersection edge, hasn't been processed
1875  // Let's start filling
1876  GridQueue *queue = new GridQueue();
1877  int total = 1;
1878 
1879  // Set to in process
1880  int mst[3];
1881  mst[0] = st[0] + vertmap[edgemap[i][0]][0] * len;
1882  mst[1] = st[1] + vertmap[edgemap[i][0]][1] * len;
1883  mst[2] = st[2] + vertmap[edgemap[i][0]][2] * len;
1884  int mdir = i / 4;
1885  setInProcessAll(mst, mdir);
1886 
1887  // Put this edge into queue
1888  queue->pushQueue(mst, mdir);
1889 
1890  // Queue processing
1891  int nst[3], dir;
1892  while (queue->popQueue(nst, dir) == 1) {
1893 #if 0
1894  dc_printf("nst: %d %d %d, dir: %d\n",
1895  nst[0] / mindimen,
1896  nst[1] / mindimen,
1897  nst[2] / mindimen,
1898  dir);
1899 #endif
1900  // locations
1901  int stMask[3][3] = {{0, 0 - len, 0 - len}, {0 - len, 0, 0 - len}, {0 - len, 0 - len, 0}};
1902  int cst[2][3];
1903  for (j = 0; j < 3; j++) {
1904  cst[0][j] = nst[j];
1905  cst[1][j] = nst[j] + stMask[dir][j];
1906  }
1907 
1908  // cells
1909  LeafNode *cs[2];
1910  for (j = 0; j < 2; j++) {
1911  cs[j] = locateLeaf(cst[j]);
1912  }
1913 
1914  // Middle sign
1915  int s = getSign(cs[0], 0);
1916 
1917  // Masks
1918  int fcCells[4] = {1, 0, 1, 0};
1919  int fcEdges[3][4][3] = {{{9, 2, 11}, {8, 1, 10}, {5, 1, 7}, {4, 2, 6}},
1920  {{10, 6, 11}, {8, 5, 9}, {1, 5, 3}, {0, 6, 2}},
1921  {{6, 10, 7}, {4, 9, 5}, {2, 9, 3}, {0, 10, 1}}};
1922 
1923  // Search for neighboring connected intersection edges
1924  for (int find = 0; find < 4; find++) {
1925  int cind = fcCells[find];
1926  int eind, edge;
1927  if (s == 0) {
1928  // Original order
1929  for (eind = 0; eind < 3; eind++) {
1930  edge = fcEdges[dir][find][eind];
1931  if (getEdgeParity(cs[cind], edge) == 1) {
1932  break;
1933  }
1934  }
1935  }
1936  else {
1937  // Inverse order
1938  for (eind = 2; eind >= 0; eind--) {
1939  edge = fcEdges[dir][find][eind];
1940  if (getEdgeParity(cs[cind], edge) == 1) {
1941  break;
1942  }
1943  }
1944  }
1945 
1946  if (eind == 3 || eind == -1) {
1947  dc_printf("Wrong! this is not a consistent sign. %d\n", eind);
1948  }
1949  else {
1950  int est[3];
1951  est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len;
1952  est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len;
1953  est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len;
1954  int edir = edge / 4;
1955 
1956  if (isInProcess(cs[cind], edge) == 0) {
1957  setInProcessAll(est, edir);
1958  queue->pushQueue(est, edir);
1959 #if 0
1960  dc_printf("Pushed: est: %d %d %d, edir: %d\n",
1961  est[0] / len,
1962  est[1] / len,
1963  est[2] / len,
1964  edir);
1965 #endif
1966  total++;
1967  }
1968  }
1969  }
1970  }
1971 
1972  dc_printf("Size of component: %d ", total);
1973 
1974  if (threshold == 0) {
1975  // Measuring stage
1976  if (total > maxtotal) {
1977  maxtotal = total;
1978  }
1979  dc_printf(".\n");
1980  delete queue;
1981  continue;
1982  }
1983 
1984  if (total >= threshold) {
1985  dc_printf("Maintained.\n");
1986  delete queue;
1987  continue;
1988  }
1989  dc_printf("Less than %d, removing...\n", threshold);
1990 
1991  // We have to remove this noise
1992 
1993  // Flip parity
1994  // setOutProcessAll(mst, mdir);
1995  flipParityAll(mst, mdir);
1996 
1997  // Put this edge into queue
1998  queue->pushQueue(mst, mdir);
1999 
2000  // Queue processing
2001  while (queue->popQueue(nst, dir) == 1) {
2002 #if 0
2003  dc_printf("nst: %d %d %d, dir: %d\n",
2004  nst[0] / mindimen,
2005  nst[1] / mindimen,
2006  nst[2] / mindimen,
2007  dir);
2008 #endif
2009  // locations
2010  int stMask[3][3] = {{0, 0 - len, 0 - len}, {0 - len, 0, 0 - len}, {0 - len, 0 - len, 0}};
2011  int cst[2][3];
2012  for (j = 0; j < 3; j++) {
2013  cst[0][j] = nst[j];
2014  cst[1][j] = nst[j] + stMask[dir][j];
2015  }
2016 
2017  // cells
2018  LeafNode *cs[2];
2019  for (j = 0; j < 2; j++)
2020  cs[j] = locateLeaf(cst[j]);
2021 
2022  // Middle sign
2023  int s = getSign(cs[0], 0);
2024 
2025  // Masks
2026  int fcCells[4] = {1, 0, 1, 0};
2027  int fcEdges[3][4][3] = {{{9, 2, 11}, {8, 1, 10}, {5, 1, 7}, {4, 2, 6}},
2028  {{10, 6, 11}, {8, 5, 9}, {1, 5, 3}, {0, 6, 2}},
2029  {{6, 10, 7}, {4, 9, 5}, {2, 9, 3}, {0, 10, 1}}};
2030 
2031  // Search for neighboring connected intersection edges
2032  for (int find = 0; find < 4; find++) {
2033  int cind = fcCells[find];
2034  int eind, edge;
2035  if (s == 0) {
2036  // Original order
2037  for (eind = 0; eind < 3; eind++) {
2038  edge = fcEdges[dir][find][eind];
2039  if (isInProcess(cs[cind], edge) == 1) {
2040  break;
2041  }
2042  }
2043  }
2044  else {
2045  // Inverse order
2046  for (eind = 2; eind >= 0; eind--) {
2047  edge = fcEdges[dir][find][eind];
2048  if (isInProcess(cs[cind], edge) == 1) {
2049  break;
2050  }
2051  }
2052  }
2053 
2054  if (eind == 3 || eind == -1) {
2055  dc_printf("Wrong! this is not a consistent sign. %d\n", eind);
2056  }
2057  else {
2058  int est[3];
2059  est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len;
2060  est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len;
2061  est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len;
2062  int edir = edge / 4;
2063 
2064  if (getEdgeParity(cs[cind], edge) == 1) {
2065  flipParityAll(est, edir);
2066  queue->pushQueue(est, edir);
2067 #if 0
2068  dc_printf("Pushed: est: %d %d %d, edir: %d\n",
2069  est[0] / len,
2070  est[1] / len,
2071  est[2] / len,
2072  edir);
2073 #endif
2074  total++;
2075  }
2076  }
2077  }
2078  }
2079 
2080  delete queue;
2081  }
2082  }
2083 
2084  return maxtotal;
2085 }
2086 
2087 int Octree::floodFill(Node *node, int st[3], int len, int height, int threshold)
2088 {
2089  int i;
2090  int maxtotal = 0;
2091 
2092  if (height == 0) {
2093  maxtotal = floodFill(&node->leaf, st, len, height, threshold);
2094  }
2095  else {
2096  // Internal cell, recur
2097  int count = 0;
2098  len >>= 1;
2099  for (i = 0; i < 8; i++) {
2100  if (node->internal.has_child(i)) {
2101  int nst[3];
2102  nst[0] = st[0] + vertmap[i][0] * len;
2103  nst[1] = st[1] + vertmap[i][1] * len;
2104  nst[2] = st[2] + vertmap[i][2] * len;
2105 
2106  int d = floodFill(node->internal.get_child(count), nst, len, height - 1, threshold);
2107  if (d > maxtotal) {
2108  maxtotal = d;
2109  }
2110  count++;
2111  }
2112  }
2113  }
2114 
2115  return maxtotal;
2116 }
2117 
2118 void Octree::writeOut()
2119 {
2120  int numQuads = 0;
2121  int numVertices = 0;
2122  int numEdges = 0;
2123 
2124  countIntersection(root, maxDepth, numQuads, numVertices, numEdges);
2125 
2126  dc_printf("Vertices counted: %d Polys counted: %d \n", numVertices, numQuads);
2127  output_mesh = alloc_output(numVertices, numQuads);
2128  int offset = 0;
2129  int st[3] = {0, 0, 0};
2130 
2131  // First, output vertices
2132  offset = 0;
2133  actualVerts = 0;
2134  actualQuads = 0;
2135 
2136  generateMinimizer(root, st, dimen, maxDepth, offset);
2137  cellProcContour(root, 0, maxDepth);
2138  dc_printf("Vertices written: %d Quads written: %d \n", offset, actualQuads);
2139 }
2140 
2141 void Octree::countIntersection(Node *node, int height, int &nedge, int &ncell, int &nface)
2142 {
2143  if (height > 0) {
2144  int total = node->internal.get_num_children();
2145  for (int i = 0; i < total; i++) {
2146  countIntersection(node->internal.get_child(i), height - 1, nedge, ncell, nface);
2147  }
2148  }
2149  else {
2150  nedge += getNumEdges2(&node->leaf);
2151 
2152  int smask = getSignMask(&node->leaf);
2153 
2154  if (use_manifold) {
2155  int comps = manifold_table[smask].comps;
2156  ncell += comps;
2157  }
2158  else {
2159  if (smask > 0 && smask < 255) {
2160  ncell++;
2161  }
2162  }
2163 
2164  for (int i = 0; i < 3; i++) {
2165  if (getFaceEdgeNum(&node->leaf, i * 2)) {
2166  nface++;
2167  }
2168  }
2169  }
2170 }
2171 
2172 /* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */
2173 static void pseudoInverse(const Eigen::Matrix3f &a, Eigen::Matrix3f &result, float tolerance)
2174 {
2175  Eigen::JacobiSVD<Eigen::Matrix3f> svd = a.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV);
2176 
2177  result = svd.matrixV() *
2178  Eigen::Vector3f((svd.singularValues().array().abs() > tolerance)
2179  .select(svd.singularValues().array().inverse(), 0))
2180  .asDiagonal() *
2181  svd.matrixU().adjoint();
2182 }
2183 
2184 static void solve_least_squares(const float halfA[],
2185  const float b[],
2186  const float midpoint[],
2187  float rvalue[])
2188 {
2189  /* calculate pseudo-inverse */
2190  Eigen::Matrix3f A, pinv;
2191  A << halfA[0], halfA[1], halfA[2], halfA[1], halfA[3], halfA[4], halfA[2], halfA[4], halfA[5];
2192  pseudoInverse(A, pinv, 0.1f);
2193 
2194  Eigen::Vector3f b2(b), mp(midpoint), result;
2195  b2 = b2 + A * -mp;
2196  result = pinv * b2 + mp;
2197 
2198  for (int i = 0; i < 3; i++)
2199  rvalue[i] = result(i);
2200 }
2201 
2202 static void mass_point(float mp[3], const float pts[12][3], const int parity[12])
2203 {
2204  int ec = 0;
2205 
2206  for (int i = 0; i < 12; i++) {
2207  if (parity[i]) {
2208  const float *p = pts[i];
2209 
2210  mp[0] += p[0];
2211  mp[1] += p[1];
2212  mp[2] += p[2];
2213 
2214  ec++;
2215  }
2216  }
2217 
2218  if (ec == 0) {
2219  return;
2220  }
2221  mp[0] /= ec;
2222  mp[1] /= ec;
2223  mp[2] /= ec;
2224 }
2225 
2226 static void minimize(float rvalue[3],
2227  float mp[3],
2228  const float pts[12][3],
2229  const float norms[12][3],
2230  const int parity[12])
2231 {
2232  float ata[6] = {0, 0, 0, 0, 0, 0};
2233  float atb[3] = {0, 0, 0};
2234  int ec = 0;
2235 
2236  for (int i = 0; i < 12; i++) {
2237  // if(getEdgeParity(leaf, i))
2238  if (parity[i]) {
2239  const float *norm = norms[i];
2240  const float *p = pts[i];
2241 
2242  // QEF
2243  ata[0] += (float)(norm[0] * norm[0]);
2244  ata[1] += (float)(norm[0] * norm[1]);
2245  ata[2] += (float)(norm[0] * norm[2]);
2246  ata[3] += (float)(norm[1] * norm[1]);
2247  ata[4] += (float)(norm[1] * norm[2]);
2248  ata[5] += (float)(norm[2] * norm[2]);
2249 
2250  const float pn = p[0] * norm[0] + p[1] * norm[1] + p[2] * norm[2];
2251 
2252  atb[0] += (float)(norm[0] * pn);
2253  atb[1] += (float)(norm[1] * pn);
2254  atb[2] += (float)(norm[2] * pn);
2255 
2256  // Minimizer
2257  mp[0] += p[0];
2258  mp[1] += p[1];
2259  mp[2] += p[2];
2260 
2261  ec++;
2262  }
2263  }
2264 
2265  if (ec == 0) {
2266  return;
2267  }
2268  mp[0] /= ec;
2269  mp[1] /= ec;
2270  mp[2] /= ec;
2271 
2272  // Solve least squares
2273  solve_least_squares(ata, atb, mp, rvalue);
2274 }
2275 
2276 void Octree::computeMinimizer(const LeafNode *leaf, int st[3], int len, float rvalue[3]) const
2277 {
2278  // First, gather all edge intersections
2279  float pts[12][3], norms[12][3];
2280  int parity[12];
2281  fillEdgeIntersections(leaf, st, len, pts, norms, parity);
2282 
2283  switch (mode) {
2284  case DUALCON_CENTROID:
2285  rvalue[0] = (float)st[0] + len / 2;
2286  rvalue[1] = (float)st[1] + len / 2;
2287  rvalue[2] = (float)st[2] + len / 2;
2288  break;
2289 
2290  case DUALCON_MASS_POINT:
2291  rvalue[0] = rvalue[1] = rvalue[2] = 0;
2292  mass_point(rvalue, pts, parity);
2293  break;
2294 
2295  default: {
2296  // Sharp features */
2297 
2298  // construct QEF and minimizer
2299  float mp[3] = {0, 0, 0};
2300  minimize(rvalue, mp, pts, norms, parity);
2301 
2302  /* Restraining the location of the minimizer */
2303  float nh1 = hermite_num * len;
2304  float nh2 = (1 + hermite_num) * len;
2305 
2306  if (rvalue[0] < st[0] - nh1 || rvalue[1] < st[1] - nh1 || rvalue[2] < st[2] - nh1 ||
2307 
2308  rvalue[0] > st[0] + nh2 || rvalue[1] > st[1] + nh2 || rvalue[2] > st[2] + nh2) {
2309  // Use mass point instead
2310  rvalue[0] = mp[0];
2311  rvalue[1] = mp[1];
2312  rvalue[2] = mp[2];
2313  }
2314  break;
2315  }
2316  }
2317 }
2318 
2319 void Octree::generateMinimizer(Node *node, int st[3], int len, int height, int &offset)
2320 {
2321  int i, j;
2322 
2323  if (height == 0) {
2324  // Leaf cell, generate
2325 
2326  // First, find minimizer
2327  float rvalue[3];
2328  rvalue[0] = (float)st[0] + len / 2;
2329  rvalue[1] = (float)st[1] + len / 2;
2330  rvalue[2] = (float)st[2] + len / 2;
2331  computeMinimizer(&node->leaf, st, len, rvalue);
2332 
2333  // Update
2334  // float fnst[3];
2335  for (j = 0; j < 3; j++) {
2336  rvalue[j] = rvalue[j] * range / dimen + origin[j];
2337  // fnst[j] = st[j] * range / dimen + origin[j];
2338  }
2339 
2340  int mult = 0, smask = getSignMask(&node->leaf);
2341 
2342  if (use_manifold) {
2343  mult = manifold_table[smask].comps;
2344  }
2345  else {
2346  if (smask > 0 && smask < 255) {
2347  mult = 1;
2348  }
2349  }
2350 
2351  for (j = 0; j < mult; j++) {
2352  add_vert(output_mesh, rvalue);
2353  }
2354 
2355  // Store the index
2356  setMinimizerIndex(&node->leaf, offset);
2357 
2358  offset += mult;
2359  }
2360  else {
2361  // Internal cell, recur
2362  int count = 0;
2363  len >>= 1;
2364  for (i = 0; i < 8; i++) {
2365  if (node->internal.has_child(i)) {
2366  int nst[3];
2367  nst[0] = st[0] + vertmap[i][0] * len;
2368  nst[1] = st[1] + vertmap[i][1] * len;
2369  nst[2] = st[2] + vertmap[i][2] * len;
2370 
2371  generateMinimizer(node->internal.get_child(count), nst, len, height - 1, offset);
2372  count++;
2373  }
2374  }
2375  }
2376 }
2377 
2378 void Octree::processEdgeWrite(Node *node[4], int /*depth*/[4], int /*maxdep*/, int dir)
2379 {
2380  // int color = 0;
2381 
2382  int i = 3;
2383  {
2384  if (getEdgeParity((LeafNode *)(node[i]), processEdgeMask[dir][i])) {
2385  int flip = 0;
2386  int edgeind = processEdgeMask[dir][i];
2387  if (getSign((LeafNode *)node[i], edgemap[edgeind][1]) > 0) {
2388  flip = 1;
2389  }
2390 
2391  int num = 0;
2392  {
2393  int ind[8];
2394  if (use_manifold) {
2395  int vind[2];
2396  int seq[4] = {0, 1, 3, 2};
2397  for (int k = 0; k < 4; k++) {
2398  getMinimizerIndices((LeafNode *)(node[seq[k]]), processEdgeMask[dir][seq[k]], vind);
2399  ind[num] = vind[0];
2400  num++;
2401 
2402  if (vind[1] != -1) {
2403  ind[num] = vind[1];
2404  num++;
2405  if (flip == 0) {
2406  ind[num - 1] = vind[0];
2407  ind[num - 2] = vind[1];
2408  }
2409  }
2410  }
2411 
2412  /* we don't use the manifold option, but if it is
2413  ever enabled again note that it can output
2414  non-quads */
2415  }
2416  else {
2417  if (flip) {
2418  ind[0] = getMinimizerIndex((LeafNode *)(node[2]));
2419  ind[1] = getMinimizerIndex((LeafNode *)(node[3]));
2420  ind[2] = getMinimizerIndex((LeafNode *)(node[1]));
2421  ind[3] = getMinimizerIndex((LeafNode *)(node[0]));
2422  }
2423  else {
2424  ind[0] = getMinimizerIndex((LeafNode *)(node[0]));
2425  ind[1] = getMinimizerIndex((LeafNode *)(node[1]));
2426  ind[2] = getMinimizerIndex((LeafNode *)(node[3]));
2427  ind[3] = getMinimizerIndex((LeafNode *)(node[2]));
2428  }
2429 
2430  add_quad(output_mesh, ind);
2431  }
2432  }
2433  return;
2434  }
2435  else {
2436  return;
2437  }
2438  }
2439 }
2440 
2441 void Octree::edgeProcContour(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir)
2442 {
2443  if (!(node[0] && node[1] && node[2] && node[3])) {
2444  return;
2445  }
2446  if (leaf[0] && leaf[1] && leaf[2] && leaf[3]) {
2447  processEdgeWrite(node, depth, maxdep, dir);
2448  }
2449  else {
2450  int i, j;
2451  Node *chd[4][8];
2452  for (j = 0; j < 4; j++) {
2453  for (i = 0; i < 8; i++) {
2454  chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2455  node[j]->internal.get_child(node[j]->internal.get_child_count(i)) :
2456  NULL;
2457  }
2458  }
2459 
2460  // 2 edge calls
2461  Node *ne[4];
2462  int le[4];
2463  int de[4];
2464  for (i = 0; i < 2; i++) {
2465  int c[4] = {edgeProcEdgeMask[dir][i][0],
2466  edgeProcEdgeMask[dir][i][1],
2467  edgeProcEdgeMask[dir][i][2],
2468  edgeProcEdgeMask[dir][i][3]};
2469 
2470  for (int j = 0; j < 4; j++) {
2471  if (leaf[j]) {
2472  le[j] = leaf[j];
2473  ne[j] = node[j];
2474  de[j] = depth[j];
2475  }
2476  else {
2477  le[j] = node[j]->internal.is_child_leaf(c[j]);
2478  ne[j] = chd[j][c[j]];
2479  de[j] = depth[j] - 1;
2480  }
2481  }
2482 
2483  edgeProcContour(ne, le, de, maxdep - 1, edgeProcEdgeMask[dir][i][4]);
2484  }
2485  }
2486 }
2487 
2488 void Octree::faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir)
2489 {
2490  if (!(node[0] && node[1])) {
2491  return;
2492  }
2493 
2494  if (!(leaf[0] && leaf[1])) {
2495  int i, j;
2496  // Fill children nodes
2497  Node *chd[2][8];
2498  for (j = 0; j < 2; j++) {
2499  for (i = 0; i < 8; i++) {
2500  chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2501  node[j]->internal.get_child(node[j]->internal.get_child_count(i)) :
2502  NULL;
2503  }
2504  }
2505 
2506  // 4 face calls
2507  Node *nf[2];
2508  int df[2];
2509  int lf[2];
2510  for (i = 0; i < 4; i++) {
2511  int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
2512  for (int j = 0; j < 2; j++) {
2513  if (leaf[j]) {
2514  lf[j] = leaf[j];
2515  nf[j] = node[j];
2516  df[j] = depth[j];
2517  }
2518  else {
2519  lf[j] = node[j]->internal.is_child_leaf(c[j]);
2520  nf[j] = chd[j][c[j]];
2521  df[j] = depth[j] - 1;
2522  }
2523  }
2524  faceProcContour(nf, lf, df, maxdep - 1, faceProcFaceMask[dir][i][2]);
2525  }
2526 
2527  // 4 edge calls
2528  int orders[2][4] = {{0, 0, 1, 1}, {0, 1, 0, 1}};
2529  Node *ne[4];
2530  int le[4];
2531  int de[4];
2532 
2533  for (i = 0; i < 4; i++) {
2534  int c[4] = {faceProcEdgeMask[dir][i][1],
2535  faceProcEdgeMask[dir][i][2],
2536  faceProcEdgeMask[dir][i][3],
2537  faceProcEdgeMask[dir][i][4]};
2538  int *order = orders[faceProcEdgeMask[dir][i][0]];
2539 
2540  for (int j = 0; j < 4; j++) {
2541  if (leaf[order[j]]) {
2542  le[j] = leaf[order[j]];
2543  ne[j] = node[order[j]];
2544  de[j] = depth[order[j]];
2545  }
2546  else {
2547  le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
2548  ne[j] = chd[order[j]][c[j]];
2549  de[j] = depth[order[j]] - 1;
2550  }
2551  }
2552 
2553  edgeProcContour(ne, le, de, maxdep - 1, faceProcEdgeMask[dir][i][5]);
2554  }
2555  }
2556 }
2557 
2558 void Octree::cellProcContour(Node *node, int leaf, int depth)
2559 {
2560  if (node == NULL) {
2561  return;
2562  }
2563 
2564  if (!leaf) {
2565  int i;
2566 
2567  // Fill children nodes
2568  Node *chd[8];
2569  for (i = 0; i < 8; i++) {
2570  chd[i] = ((!leaf) && node->internal.has_child(i)) ?
2571  node->internal.get_child(node->internal.get_child_count(i)) :
2572  NULL;
2573  }
2574 
2575  // 8 Cell calls
2576  for (i = 0; i < 8; i++) {
2577  cellProcContour(chd[i], node->internal.is_child_leaf(i), depth - 1);
2578  }
2579 
2580  // 12 face calls
2581  Node *nf[2];
2582  int lf[2];
2583  int df[2] = {depth - 1, depth - 1};
2584  for (i = 0; i < 12; i++) {
2585  int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
2586 
2587  lf[0] = node->internal.is_child_leaf(c[0]);
2588  lf[1] = node->internal.is_child_leaf(c[1]);
2589 
2590  nf[0] = chd[c[0]];
2591  nf[1] = chd[c[1]];
2592 
2593  faceProcContour(nf, lf, df, depth - 1, cellProcFaceMask[i][2]);
2594  }
2595 
2596  // 6 edge calls
2597  Node *ne[4];
2598  int le[4];
2599  int de[4] = {depth - 1, depth - 1, depth - 1, depth - 1};
2600  for (i = 0; i < 6; i++) {
2601  int c[4] = {cellProcEdgeMask[i][0],
2602  cellProcEdgeMask[i][1],
2603  cellProcEdgeMask[i][2],
2604  cellProcEdgeMask[i][3]};
2605 
2606  for (int j = 0; j < 4; j++) {
2607  le[j] = node->internal.is_child_leaf(c[j]);
2608  ne[j] = chd[c[j]];
2609  }
2610 
2611  edgeProcContour(ne, le, de, depth - 1, cellProcEdgeMask[i][4]);
2612  }
2613  }
2614 }
2615 
2616 void Octree::processEdgeParity(LeafNode *node[4], int /*depth*/[4], int /*maxdep*/, int dir)
2617 {
2618  int con = 0;
2619  for (int i = 0; i < 4; i++) {
2620  // Minimal cell
2621  // if(op == 0)
2622  {
2623  if (getEdgeParity(node[i], processEdgeMask[dir][i])) {
2624  con = 1;
2625  break;
2626  }
2627  }
2628  }
2629 
2630  if (con == 1) {
2631  for (int i = 0; i < 4; i++) {
2632  setEdge(node[i], processEdgeMask[dir][i]);
2633  }
2634  }
2635 }
2636 
2637 void Octree::edgeProcParity(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir)
2638 {
2639  if (!(node[0] && node[1] && node[2] && node[3])) {
2640  return;
2641  }
2642  if (leaf[0] && leaf[1] && leaf[2] && leaf[3]) {
2643  processEdgeParity((LeafNode **)node, depth, maxdep, dir);
2644  }
2645  else {
2646  int i, j;
2647  Node *chd[4][8];
2648  for (j = 0; j < 4; j++) {
2649  for (i = 0; i < 8; i++) {
2650  chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2651  node[j]->internal.get_child(node[j]->internal.get_child_count(i)) :
2652  NULL;
2653  }
2654  }
2655 
2656  // 2 edge calls
2657  Node *ne[4];
2658  int le[4];
2659  int de[4];
2660  for (i = 0; i < 2; i++) {
2661  int c[4] = {edgeProcEdgeMask[dir][i][0],
2662  edgeProcEdgeMask[dir][i][1],
2663  edgeProcEdgeMask[dir][i][2],
2664  edgeProcEdgeMask[dir][i][3]};
2665 
2666  // int allleaf = 1;
2667  for (int j = 0; j < 4; j++) {
2668 
2669  if (leaf[j]) {
2670  le[j] = leaf[j];
2671  ne[j] = node[j];
2672  de[j] = depth[j];
2673  }
2674  else {
2675  le[j] = node[j]->internal.is_child_leaf(c[j]);
2676  ne[j] = chd[j][c[j]];
2677  de[j] = depth[j] - 1;
2678  }
2679  }
2680 
2681  edgeProcParity(ne, le, de, maxdep - 1, edgeProcEdgeMask[dir][i][4]);
2682  }
2683  }
2684 }
2685 
2686 void Octree::faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir)
2687 {
2688  if (!(node[0] && node[1])) {
2689  return;
2690  }
2691 
2692  if (!(leaf[0] && leaf[1])) {
2693  int i, j;
2694  // Fill children nodes
2695  Node *chd[2][8];
2696  for (j = 0; j < 2; j++) {
2697  for (i = 0; i < 8; i++) {
2698  chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2699  node[j]->internal.get_child(node[j]->internal.get_child_count(i)) :
2700  NULL;
2701  }
2702  }
2703 
2704  // 4 face calls
2705  Node *nf[2];
2706  int df[2];
2707  int lf[2];
2708  for (i = 0; i < 4; i++) {
2709  int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
2710  for (int j = 0; j < 2; j++) {
2711  if (leaf[j]) {
2712  lf[j] = leaf[j];
2713  nf[j] = node[j];
2714  df[j] = depth[j];
2715  }
2716  else {
2717  lf[j] = node[j]->internal.is_child_leaf(c[j]);
2718  nf[j] = chd[j][c[j]];
2719  df[j] = depth[j] - 1;
2720  }
2721  }
2722  faceProcParity(nf, lf, df, maxdep - 1, faceProcFaceMask[dir][i][2]);
2723  }
2724 
2725  // 4 edge calls
2726  int orders[2][4] = {{0, 0, 1, 1}, {0, 1, 0, 1}};
2727  Node *ne[4];
2728  int le[4];
2729  int de[4];
2730 
2731  for (i = 0; i < 4; i++) {
2732  int c[4] = {faceProcEdgeMask[dir][i][1],
2733  faceProcEdgeMask[dir][i][2],
2734  faceProcEdgeMask[dir][i][3],
2735  faceProcEdgeMask[dir][i][4]};
2736  int *order = orders[faceProcEdgeMask[dir][i][0]];
2737 
2738  for (int j = 0; j < 4; j++) {
2739  if (leaf[order[j]]) {
2740  le[j] = leaf[order[j]];
2741  ne[j] = node[order[j]];
2742  de[j] = depth[order[j]];
2743  }
2744  else {
2745  le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
2746  ne[j] = chd[order[j]][c[j]];
2747  de[j] = depth[order[j]] - 1;
2748  }
2749  }
2750 
2751  edgeProcParity(ne, le, de, maxdep - 1, faceProcEdgeMask[dir][i][5]);
2752  }
2753  }
2754 }
2755 
2756 void Octree::cellProcParity(Node *node, int leaf, int depth)
2757 {
2758  if (node == NULL) {
2759  return;
2760  }
2761 
2762  if (!leaf) {
2763  int i;
2764 
2765  // Fill children nodes
2766  Node *chd[8];
2767  for (i = 0; i < 8; i++) {
2768  chd[i] = ((!leaf) && node->internal.has_child(i)) ?
2769  node->internal.get_child(node->internal.get_child_count(i)) :
2770  NULL;
2771  }
2772 
2773  // 8 Cell calls
2774  for (i = 0; i < 8; i++) {
2775  cellProcParity(chd[i], node->internal.is_child_leaf(i), depth - 1);
2776  }
2777 
2778  // 12 face calls
2779  Node *nf[2];
2780  int lf[2];
2781  int df[2] = {depth - 1, depth - 1};
2782  for (i = 0; i < 12; i++) {
2783  int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
2784 
2785  lf[0] = node->internal.is_child_leaf(c[0]);
2786  lf[1] = node->internal.is_child_leaf(c[1]);
2787 
2788  nf[0] = chd[c[0]];
2789  nf[1] = chd[c[1]];
2790 
2791  faceProcParity(nf, lf, df, depth - 1, cellProcFaceMask[i][2]);
2792  }
2793 
2794  // 6 edge calls
2795  Node *ne[4];
2796  int le[4];
2797  int de[4] = {depth - 1, depth - 1, depth - 1, depth - 1};
2798  for (i = 0; i < 6; i++) {
2799  int c[4] = {cellProcEdgeMask[i][0],
2800  cellProcEdgeMask[i][1],
2801  cellProcEdgeMask[i][2],
2802  cellProcEdgeMask[i][3]};
2803 
2804  for (int j = 0; j < 4; j++) {
2805  le[j] = node->internal.is_child_leaf(c[j]);
2806  ne[j] = chd[c[j]];
2807  }
2808 
2809  edgeProcParity(ne, le, de, depth - 1, cellProcEdgeMask[i][4]);
2810  }
2811  }
2812 }
2813 
2814 /* definitions for global arrays */
2815 const int edgemask[3] = {5, 3, 6};
2816 
2817 const int faceMap[6][4] = {
2818  {4, 8, 5, 9},
2819  {6, 10, 7, 11},
2820  {0, 8, 1, 10},
2821  {2, 9, 3, 11},
2822  {0, 4, 2, 6},
2823  {1, 5, 3, 7},
2824 };
2825 
2826 const int cellProcFaceMask[12][3] = {
2827  {0, 4, 0},
2828  {1, 5, 0},
2829  {2, 6, 0},
2830  {3, 7, 0},
2831  {0, 2, 1},
2832  {4, 6, 1},
2833  {1, 3, 1},
2834  {5, 7, 1},
2835  {0, 1, 2},
2836  {2, 3, 2},
2837  {4, 5, 2},
2838  {6, 7, 2},
2839 };
2840 
2841 const int cellProcEdgeMask[6][5] = {
2842  {0, 1, 2, 3, 0},
2843  {4, 5, 6, 7, 0},
2844  {0, 4, 1, 5, 1},
2845  {2, 6, 3, 7, 1},
2846  {0, 2, 4, 6, 2},
2847  {1, 3, 5, 7, 2},
2848 };
2849 
2850 const int faceProcFaceMask[3][4][3] = {{{4, 0, 0}, {5, 1, 0}, {6, 2, 0}, {7, 3, 0}},
2851  {{2, 0, 1}, {6, 4, 1}, {3, 1, 1}, {7, 5, 1}},
2852  {{1, 0, 2}, {3, 2, 2}, {5, 4, 2}, {7, 6, 2}}};
2853 
2854 const int faceProcEdgeMask[3][4][6] = {
2855  {{1, 4, 0, 5, 1, 1}, {1, 6, 2, 7, 3, 1}, {0, 4, 6, 0, 2, 2}, {0, 5, 7, 1, 3, 2}},
2856  {{0, 2, 3, 0, 1, 0}, {0, 6, 7, 4, 5, 0}, {1, 2, 0, 6, 4, 2}, {1, 3, 1, 7, 5, 2}},
2857  {{1, 1, 0, 3, 2, 0}, {1, 5, 4, 7, 6, 0}, {0, 1, 5, 0, 4, 1}, {0, 3, 7, 2, 6, 1}}};
2858 
2859 const int edgeProcEdgeMask[3][2][5] = {
2860  {{3, 2, 1, 0, 0}, {7, 6, 5, 4, 0}},
2861  {{5, 1, 4, 0, 1}, {7, 3, 6, 2, 1}},
2862  {{6, 4, 2, 0, 2}, {7, 5, 3, 1, 2}},
2863 };
2864 
2865 const int processEdgeMask[3][4] = {
2866  {3, 2, 1, 0},
2867  {7, 5, 6, 4},
2868  {11, 10, 9, 8},
2869 };
2870 
2871 const int dirCell[3][4][3] = {{{0, -1, -1}, {0, -1, 0}, {0, 0, -1}, {0, 0, 0}},
2872  {{-1, 0, -1}, {-1, 0, 0}, {0, 0, -1}, {0, 0, 0}},
2873  {{-1, -1, 0}, {-1, 0, 0}, {0, -1, 0}, {0, 0, 0}}};
2874 
2875 const int dirEdge[3][4] = {
2876  {3, 2, 1, 0},
2877  {7, 6, 5, 4},
2878  {11, 10, 9, 8},
2879 };
2880 
typedef float(TangentPoint)[2]
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble y1
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei height
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint y
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble GLdouble x2
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint order
Read Guarded memory(de)allocation.
const int edgemap[12][2]
Definition: Projections.cpp:28
const int vertmap[8][3]
Definition: Projections.cpp:10
#define GRID_DIMENSION
Definition: Projections.h:10
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
#define A
static DBVT_INLINE btScalar size(const btDbvtVolume &a)
Definition: btDbvt.cpp:52
int numVertices() const
SIMD_FORCE_INLINE void mult(const btTransform &t1, const btTransform &t2)
Set the current transform as the value of the product of two transforms.
Definition: btTransform.h:76
SIMD_FORCE_INLINE btScalar norm() const
Return the norm (length) of the vector.
Definition: btVector3.h:263
void shift(int off[3])
int isIntersecting() const
int isIntersectingPrimary(int edgeInd) const
unsigned char getBoxMask()
TriangleProjection * inherit
Inheritable portion.
Definition: Projections.h:72
float getIntersectionPrimary(int edgeInd) const
Definition: cubes.h:9
virtual float getBoundingBox(float origin[3])=0
Get bounding box.
virtual Triangle * getNextTriangle()=0
Get next triangle.
virtual int getNumTriangles()=0
Get number of triangles.
int use_flood_fill
Definition: octree.h:247
int use_manifold
Definition: octree.h:250
int dimen
Definition: octree.h:224
float range
Definition: octree.h:232
int maxTrianglePerCell
Definition: octree.h:243
float thresh
Definition: octree.h:248
ModelReader * reader
Definition: octree.h:218
VirtualMemoryAllocator * alloc[9]
Definition: octree.h:211
int maxDepth
Definition: octree.h:228
int nodeSpace
Definition: octree.h:236
VirtualMemoryAllocator * leafalloc[4]
Definition: octree.h:212
~Octree()
Definition: octree.cpp:91
float origin[3]
Definition: octree.h:231
Cubes * cubes
Definition: octree.h:221
int nodeCount
Definition: octree.h:235
float hermite_num
Definition: octree.h:252
int minshift
Definition: octree.h:225
int actualQuads
Definition: octree.h:239
DualConMode mode
Definition: octree.h:254
void scanConvert()
Definition: octree.cpp:97
int mindimen
Definition: octree.h:225
Octree(ModelReader *mr, DualConAllocOutput alloc_output_func, DualConAddVert add_vert_func, DualConAddQuad add_quad_func, DualConFlags flags, DualConMode mode, int depth, float threshold, float hermite_num)
Definition: octree.cpp:31
Node * root
Definition: octree.h:215
int actualVerts
Definition: octree.h:239
virtual int getBytes()=0
virtual int getAllocated()=0
virtual int getAll()=0
virtual void destroy()=0
virtual void printInfo()=0
OperationNode * node
SyclQueue * queue
SyclQueue void void size_t num_bytes void
int len
Definition: draw_manager.c:108
DualConMode
Definition: dualcon.h:47
@ DUALCON_CENTROID
Definition: dualcon.h:49
@ DUALCON_MASS_POINT
Definition: dualcon.h:51
void(* DualConAddQuad)(void *output, const int vert_indices[4])
Definition: dualcon.h:41
void(* DualConAddVert)(void *output, const float co[3])
Definition: dualcon.h:39
DualConFlags
Definition: dualcon.h:43
@ DUALCON_FLOOD_FILL
Definition: dualcon.h:44
void *(* DualConAllocOutput)(int totvert, int totquad)
Definition: dualcon.h:37
uint pos
int count
ccl_gpu_kernel_postfix ccl_global float int int int int float threshold
ccl_gpu_kernel_postfix ccl_global float int int sy
ccl_gpu_kernel_postfix ccl_global float int int int int float bool int offset
ccl_gpu_kernel_postfix ccl_global float int sx
const ManifoldIndices manifold_table[256]
ccl_device_inline float4 mask(const int4 &mask, const float4 &a)
Definition: math_float4.h:513
static ulong * next
static unsigned c
Definition: RandGen.cpp:83
static unsigned a[3]
Definition: RandGen.cpp:78
T midpoint(const T &a, const T &b)
SymEdge< T > * prev(const SymEdge< T > *se)
Definition: delaunay_2d.cc:105
static const pxr::TfToken st("st", pxr::TfToken::Immortal)
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
const int processEdgeMask[3][4]
Definition: octree.cpp:2865
const int edgeProcEdgeMask[3][2][5]
Definition: octree.cpp:2859
static void pseudoInverse(const Eigen::Matrix3f &a, Eigen::Matrix3f &result, float tolerance)
Definition: octree.cpp:2173
const int faceProcFaceMask[3][4][3]
Definition: octree.cpp:2850
static void minimize(float rvalue[3], float mp[3], const float pts[12][3], const float norms[12][3], const int parity[12])
Definition: octree.cpp:2226
static void solve_least_squares(const float halfA[], const float b[], const float midpoint[], float rvalue[])
Definition: octree.cpp:2184
const int cellProcFaceMask[12][3]
Definition: octree.cpp:2826
const int faceMap[6][4]
Definition: octree.cpp:2817
const int faceProcEdgeMask[3][4][6]
Definition: octree.cpp:2854
#define dc_printf(...)
Definition: octree.cpp:26
const int cellProcEdgeMask[6][5]
Definition: octree.cpp:2841
const int dirEdge[3][4]
Definition: octree.cpp:2875
static void mass_point(float mp[3], const float pts[12][3], const int parity[12])
Definition: octree.cpp:2202
const int dirCell[3][4][3]
Definition: octree.cpp:2871
const int edgemask[3]
Definition: octree.cpp:2815
__int64 int64_t
Definition: stdint.h:89
Node * get_child(int count)
Definition: octree.h:67
static int childrenCountTable[256][8]
Definition: octree.h:43
static int numChildrenTable[256]
Definition: octree.h:42
int has_child(int index) const
Definition: octree.h:61
static int childrenIndexTable[256][8]
Definition: octree.h:44
void fill_children(Node *children[8], int leaf[8])
Definition: octree.h:98
int is_child_leaf(int index) const
Definition: octree.h:55
LeafNode leaf
Definition: octree.h:164
InternalNode internal
Definition: octree.h:163
PathElement * next
Definition: octree.h:188
int pos[3]
Definition: octree.h:185
int length
Definition: octree.h:197
PathElement * tail
Definition: octree.h:194
PathList * next
Definition: octree.h:200
PathElement * head
Definition: octree.h:193
double norm[3]
Normal of the triangle.
Definition: Projections.h:52
float vt[3][3]
Definition: GeoCommon.h:29
SymEdge< Arith_t > * next
Definition: delaunay_2d.cc:83