dd708b27019faeff9fdf47a88443138e209f984a
[blender-staging.git] / source / blender / blenlib / intern / graph.c
1 /**
2  * $Id:
3  *
4  * ***** BEGIN GPL LICENSE BLOCK *****
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software Foundation,
18  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19  *
20  * Contributor(s): Martin Poirier
21  *
22  * ***** END GPL LICENSE BLOCK *****
23  * graph.c: Common graph interface and methods
24  */
25
26 #include <float.h> 
27 #include <math.h>
28
29 #include "MEM_guardedalloc.h"
30
31 #include "BLI_graph.h"
32 #include "BLI_blenlib.h"
33 #include "BLI_math.h"
34
35 #include "BKE_utildefines.h"
36
37 static void testRadialSymmetry(BGraph *graph, BNode* root_node, RadialArc* ring, int total, float axis[3], float limit, int group);
38
39 static void handleAxialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit);
40 static void testAxialSymmetry(BGraph *graph, BNode* root_node, BNode* node1, BNode* node2, BArc* arc1, BArc* arc2, float axis[3], float limit, int group);
41 static void flagAxialSymmetry(BNode *root_node, BNode *end_node, BArc *arc, int group);
42
43 void BLI_freeNode(BGraph *graph, BNode *node)
44 {
45         if (node->arcs)
46         {
47                 MEM_freeN(node->arcs);
48         }
49         
50         if (graph->free_node)
51         {
52                 graph->free_node(node);
53         }
54 }
55
56 void BLI_removeNode(BGraph *graph, BNode *node)
57 {
58         BLI_freeNode(graph, node);
59         BLI_freelinkN(&graph->nodes, node);
60 }
61
62 BNode *BLI_otherNode(BArc *arc, BNode *node)
63 {
64         return (arc->head == node) ? arc->tail : arc->head;
65 }
66
67 void BLI_removeArc(BGraph *graph, BArc *arc)
68 {
69         if (graph->free_arc)
70         {
71                 graph->free_arc(arc);
72         }
73
74         BLI_freelinkN(&graph->arcs, arc);
75 }
76
77 void BLI_flagNodes(BGraph *graph, int flag)
78 {
79         BNode *node;
80         
81         for(node = graph->nodes.first; node; node = node->next)
82         {
83                 node->flag = flag;
84         }
85 }
86
87 void BLI_flagArcs(BGraph *graph, int flag)
88 {
89         BArc *arc;
90         
91         for(arc = graph->arcs.first; arc; arc = arc->next)
92         {
93                 arc->flag = flag;
94         }
95 }
96
97 static void addArcToNodeAdjacencyList(BNode *node, BArc *arc)
98 {
99         node->arcs[node->flag] = arc;
100         node->flag++;
101 }
102
103 void BLI_buildAdjacencyList(BGraph *graph)
104 {
105         BNode *node;
106         BArc *arc;
107
108         for(node = graph->nodes.first; node; node = node->next)
109         {
110                 if (node->arcs != NULL)
111                 {
112                         MEM_freeN(node->arcs);
113                 }
114                 
115                 node->arcs = MEM_callocN((node->degree) * sizeof(BArc*), "adjacency list");
116                 
117                 /* temporary use to indicate the first index available in the lists */
118                 node->flag = 0;
119         }
120
121         for(arc = graph->arcs.first; arc; arc= arc->next)
122         {
123                 addArcToNodeAdjacencyList(arc->head, arc);
124                 addArcToNodeAdjacencyList(arc->tail, arc);
125         }
126
127         for(node = graph->nodes.first; node; node = node->next)
128         {
129                 if (node->degree != node->flag)
130                 {
131                         printf("error in node [%p]. Added only %i arcs out of %i\n", node, node->flag, node->degree);
132                 }
133         }
134 }
135
136 void BLI_rebuildAdjacencyListForNode(BGraph* graph, BNode *node)
137 {
138         BArc *arc;
139
140         if (node->arcs != NULL)
141         {
142                 MEM_freeN(node->arcs);
143         }
144         
145         node->arcs = MEM_callocN((node->degree) * sizeof(BArc*), "adjacency list");
146         
147         /* temporary use to indicate the first index available in the lists */
148         node->flag = 0;
149
150         for(arc = graph->arcs.first; arc; arc= arc->next)
151         {
152                 if (arc->head == node)
153                 {
154                         addArcToNodeAdjacencyList(arc->head, arc);
155                 }
156                 else if (arc->tail == node)
157                 {
158                         addArcToNodeAdjacencyList(arc->tail, arc);
159                 }
160         }
161
162         if (node->degree != node->flag)
163         {
164                 printf("error in node [%p]. Added only %i arcs out of %i\n", node, node->flag, node->degree);
165         }
166 }
167
168 void BLI_freeAdjacencyList(BGraph *graph)
169 {
170         BNode *node;
171
172         for(node = graph->nodes.first; node; node = node->next)
173         {
174                 if (node->arcs != NULL)
175                 {
176                         MEM_freeN(node->arcs);
177                         node->arcs = NULL;
178                 }
179         }
180 }
181
182 int BLI_hasAdjacencyList(BGraph *graph)
183 {
184         BNode *node;
185         
186         for(node = graph->nodes.first; node; node = node->next)
187         {
188                 if (node->arcs == NULL)
189                 {
190                         return 0;
191                 }
192         }
193         
194         return 1;
195 }
196
197 void BLI_replaceNodeInArc(BGraph *graph, BArc *arc, BNode *node_src, BNode *node_replaced)
198 {
199         if (arc->head == node_replaced)
200         {
201                 arc->head = node_src;
202                 node_src->degree++;
203         }
204
205         if (arc->tail == node_replaced)
206         {
207                 arc->tail = node_src;
208                 node_src->degree++;
209         }
210         
211         if (arc->head == arc->tail)
212         {
213                 node_src->degree -= 2;
214                 
215                 graph->free_arc(arc);
216                 BLI_freelinkN(&graph->arcs, arc);
217         }
218
219         if (node_replaced->degree == 0)
220         {
221                 BLI_removeNode(graph, node_replaced);
222         }
223 }
224
225 void BLI_replaceNode(BGraph *graph, BNode *node_src, BNode *node_replaced)
226 {
227         BArc *arc, *next_arc;
228         
229         for (arc = graph->arcs.first; arc; arc = next_arc)
230         {
231                 next_arc = arc->next;
232                 
233                 if (arc->head == node_replaced)
234                 {
235                         arc->head = node_src;
236                         node_replaced->degree--;
237                         node_src->degree++;
238                 }
239
240                 if (arc->tail == node_replaced)
241                 {
242                         arc->tail = node_src;
243                         node_replaced->degree--;
244                         node_src->degree++;
245                 }
246                 
247                 if (arc->head == arc->tail)
248                 {
249                         node_src->degree -= 2;
250                         
251                         graph->free_arc(arc);
252                         BLI_freelinkN(&graph->arcs, arc);
253                 }
254         }
255         
256         if (node_replaced->degree == 0)
257         {
258                 BLI_removeNode(graph, node_replaced);
259         }
260 }
261
262 void BLI_removeDoubleNodes(BGraph *graph, float limit)
263 {
264         BNode *node_src, *node_replaced;
265         
266         for(node_src = graph->nodes.first; node_src; node_src = node_src->next)
267         {
268                 for(node_replaced = graph->nodes.first; node_replaced; node_replaced = node_replaced->next)
269                 {
270                         if (node_replaced != node_src && len_v3v3(node_replaced->p, node_src->p) <= limit)
271                         {
272                                 BLI_replaceNode(graph, node_src, node_replaced);
273                         }
274                 }
275         }
276         
277 }
278
279 BNode * BLI_FindNodeByPosition(BGraph *graph, float *p, float limit)
280 {
281         BNode *closest_node = NULL, *node;
282         float min_distance = 0.0f;
283         
284         for(node = graph->nodes.first; node; node = node->next)
285         {
286                 float distance = len_v3v3(p, node->p); 
287                 if (distance <= limit && (closest_node == NULL || distance < min_distance))
288                 {
289                         closest_node = node;
290                         min_distance = distance;
291                 }
292         }
293         
294         return closest_node;
295 }
296 /************************************* SUBGRAPH DETECTION **********************************************/
297
298 void flagSubgraph(BNode *node, int subgraph)
299 {
300         if (node->subgraph_index == 0)
301         {
302                 BArc *arc;
303                 int i;
304                 
305                 node->subgraph_index = subgraph;
306                 
307                 for(i = 0; i < node->degree; i++)
308                 {
309                         arc = node->arcs[i];
310                         flagSubgraph(BLI_otherNode(arc, node), subgraph);
311                 }
312         }
313
314
315 int BLI_FlagSubgraphs(BGraph *graph)
316 {
317         BNode *node;
318         int subgraph = 0;
319
320         if (BLI_hasAdjacencyList(graph) == 0)
321         {
322                 BLI_buildAdjacencyList(graph);
323         }
324         
325         for(node = graph->nodes.first; node; node = node->next)
326         {
327                 node->subgraph_index = 0;
328         }
329         
330         for (node = graph->nodes.first; node; node = node->next)
331         {
332                 if (node->subgraph_index == 0)
333                 {
334                         subgraph++;
335                         flagSubgraph(node, subgraph);
336                 }
337         }
338         
339         return subgraph;
340 }
341
342 void BLI_ReflagSubgraph(BGraph *graph, int old_subgraph, int new_subgraph)
343 {
344         BNode *node;
345
346         for (node = graph->nodes.first; node; node = node->next)
347         {
348                 if (node->flag == old_subgraph)
349                 {
350                         node->flag = new_subgraph;
351                 }
352         }
353 }
354
355 /*************************************** CYCLE DETECTION ***********************************************/
356
357 static int detectCycle(BNode *node, BArc *src_arc)
358 {
359         int value = 0;
360         
361         if (node->flag == 0)
362         {
363                 int i;
364
365                 /* mark node as visited */
366                 node->flag = 1;
367
368                 for(i = 0; i < node->degree && value == 0; i++)
369                 {
370                         BArc *arc = node->arcs[i];
371                         
372                         /* don't go back on the source arc */
373                         if (arc != src_arc)
374                         {
375                                 value = detectCycle(BLI_otherNode(arc, node), arc);
376                         }
377                 }
378         }
379         else
380         {
381                 value = 1;
382         }
383         
384         return value;
385 }
386
387 int     BLI_isGraphCyclic(BGraph *graph)
388 {
389         BNode *node;
390         int value = 0;
391         
392         /* NEED TO CHECK IF ADJACENCY LIST EXIST */
393         
394         /* Mark all nodes as not visited */
395         BLI_flagNodes(graph, 0);
396
397         /* detectCycles in subgraphs */ 
398         for(node = graph->nodes.first; node && value == 0; node = node->next)
399         {
400                 /* only for nodes in subgraphs that haven't been visited yet */
401                 if (node->flag == 0)
402                 {
403                         value = value || detectCycle(node, NULL);
404                 }               
405         }
406         
407         return value;
408 }
409
410 BArc * BLI_findConnectedArc(BGraph *graph, BArc *arc, BNode *v)
411 {
412         BArc *nextArc = arc->next;
413         
414         for(nextArc = graph->arcs.first; nextArc; nextArc = nextArc->next)
415         {
416                 if (arc != nextArc && (nextArc->head == v || nextArc->tail == v))
417                 {
418                         break;
419                 }
420         }
421         
422         return nextArc;
423 }
424
425 /*********************************** GRAPH AS TREE FUNCTIONS *******************************************/
426
427 int subtreeShape(BNode *node, BArc *rootArc, int include_root)
428 {
429         int depth = 0;
430         
431         node->flag = 1;
432         
433         if (include_root)
434         {
435                 BNode *newNode = BLI_otherNode(rootArc, node);
436                 return subtreeShape(newNode, rootArc, 0);
437         }
438         else
439         {
440                 /* Base case, no arcs leading away */
441                 if (node->arcs == NULL || *(node->arcs) == NULL)
442                 {
443                         return 0;
444                 }
445                 else
446                 {
447                         int i;
448         
449                         for(i = 0; i < node->degree; i++)
450                         {
451                                 BArc *arc = node->arcs[i];
452                                 BNode *newNode = BLI_otherNode(arc, node);
453                                 
454                                 /* stop immediate and cyclic backtracking */
455                                 if (arc != rootArc && newNode->flag == 0)
456                                 {
457                                         depth += subtreeShape(newNode, arc, 0);
458                                 }
459                         }
460                 }
461                 
462                 return SHAPE_RADIX * depth + 1;
463         }
464 }
465
466 int BLI_subtreeShape(BGraph *graph, BNode *node, BArc *rootArc, int include_root)
467 {
468         BLI_flagNodes(graph, 0);
469         return subtreeShape(node, rootArc, include_root);
470 }
471
472 float BLI_subtreeLength(BNode *node)
473 {
474         float length = 0;
475         int i;
476
477         node->flag = 0; /* flag node as visited */
478
479         for(i = 0; i < node->degree; i++)
480         {
481                 BArc *arc = node->arcs[i];
482                 BNode *other_node = BLI_otherNode(arc, node);
483                 
484                 if (other_node->flag != 0)
485                 {
486                         float subgraph_length = arc->length + BLI_subtreeLength(other_node); 
487                         length = MAX2(length, subgraph_length);
488                 }
489         }
490         
491         return length;
492 }
493
494 void BLI_calcGraphLength(BGraph *graph)
495 {
496         float length = 0;
497         int nb_subgraphs;
498         int i;
499         
500         nb_subgraphs = BLI_FlagSubgraphs(graph);
501         
502         for (i = 1; i <= nb_subgraphs; i++)
503         {
504                 BNode *node;
505                 
506                 for (node = graph->nodes.first; node; node = node->next)
507                 {
508                         /* start on an external node  of the subgraph */
509                         if (node->subgraph_index == i && node->degree == 1)
510                         {
511                                 float subgraph_length = BLI_subtreeLength(node);
512                                 length = MAX2(length, subgraph_length);
513                                 break;
514                         }
515                 }
516         }
517         
518         graph->length = length;
519 }
520
521 /********************************* SYMMETRY DETECTION **************************************************/
522
523 static void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit);
524
525 void BLI_mirrorAlongAxis(float v[3], float center[3], float axis[3])
526 {
527         float dv[3], pv[3];
528         
529         sub_v3_v3v3(dv, v, center);
530         project_v3_v3v3(pv, dv, axis);
531         mul_v3_fl(pv, -2);
532         add_v3_v3v3(v, v, pv);
533 }
534
535 static void testRadialSymmetry(BGraph *graph, BNode* root_node, RadialArc* ring, int total, float axis[3], float limit, int group)
536 {
537         int symmetric = 1;
538         int i;
539         
540         /* sort ring by angle */
541         for (i = 0; i < total - 1; i++)
542         {
543                 float minAngle = FLT_MAX;
544                 int minIndex = -1;
545                 int j;
546
547                 for (j = i + 1; j < total; j++)
548                 {
549                         float angle = dot_v3v3(ring[i].n, ring[j].n);
550
551                         /* map negative values to 1..2 */
552                         if (angle < 0)
553                         {
554                                 angle = 1 - angle;
555                         }
556
557                         if (angle < minAngle)
558                         {
559                                 minIndex = j;
560                                 minAngle = angle;
561                         }
562                 }
563
564                 /* swap if needed */
565                 if (minIndex != i + 1)
566                 {
567                         RadialArc tmp;
568                         tmp = ring[i + 1];
569                         ring[i + 1] = ring[minIndex];
570                         ring[minIndex] = tmp;
571                 }
572         }
573
574         for (i = 0; i < total && symmetric; i++)
575         {
576                 BNode *node1, *node2;
577                 float tangent[3];
578                 float normal[3];
579                 float p[3];
580                 int j = (i + 1) % total; /* next arc in the circular list */
581
582                 add_v3_v3v3(tangent, ring[i].n, ring[j].n);
583                 cross_v3_v3v3(normal, tangent, axis);
584                 
585                 node1 = BLI_otherNode(ring[i].arc, root_node);
586                 node2 = BLI_otherNode(ring[j].arc, root_node);
587
588                 VECCOPY(p, node2->p);
589                 BLI_mirrorAlongAxis(p, root_node->p, normal);
590                 
591                 /* check if it's within limit before continuing */
592                 if (len_v3v3(node1->p, p) > limit)
593                 {
594                         symmetric = 0;
595                 }
596
597         }
598
599         if (symmetric)
600         {
601                 /* mark node as symmetric physically */
602                 VECCOPY(root_node->symmetry_axis, axis);
603                 root_node->symmetry_flag |= SYM_PHYSICAL;
604                 root_node->symmetry_flag |= SYM_RADIAL;
605                 
606                 /* FLAG SYMMETRY GROUP */
607                 for (i = 0; i < total; i++)
608                 {
609                         ring[i].arc->symmetry_group = group;
610                         ring[i].arc->symmetry_flag = SYM_SIDE_RADIAL + i;
611                 }
612
613                 if (graph->radial_symmetry)
614                 {
615                         graph->radial_symmetry(root_node, ring, total);
616                 }
617         }
618 }
619
620 static void handleRadialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit)
621 {
622         RadialArc *ring = NULL;
623         RadialArc *unit;
624         int total = 0;
625         int group;
626         int first;
627         int i;
628
629         /* mark topological symmetry */
630         root_node->symmetry_flag |= SYM_TOPOLOGICAL;
631
632         /* total the number of arcs in the symmetry ring */
633         for (i = 0; i < root_node->degree; i++)
634         {
635                 BArc *connectedArc = root_node->arcs[i];
636                 
637                 /* depth is store as a negative in flag. symmetry level is positive */
638                 if (connectedArc->symmetry_level == -depth)
639                 {
640                         total++;
641                 }
642         }
643
644         ring = MEM_callocN(sizeof(RadialArc) * total, "radial symmetry ring");
645         unit = ring;
646
647         /* fill in the ring */
648         for (unit = ring, i = 0; i < root_node->degree; i++)
649         {
650                 BArc *connectedArc = root_node->arcs[i];
651                 
652                 /* depth is store as a negative in flag. symmetry level is positive */
653                 if (connectedArc->symmetry_level == -depth)
654                 {
655                         BNode *otherNode = BLI_otherNode(connectedArc, root_node);
656                         float vec[3];
657
658                         unit->arc = connectedArc;
659
660                         /* project the node to node vector on the symmetry plane */
661                         sub_v3_v3v3(unit->n, otherNode->p, root_node->p);
662                         project_v3_v3v3(vec, unit->n, axis);
663                         sub_v3_v3v3(unit->n, unit->n, vec);
664
665                         normalize_v3(unit->n);
666
667                         unit++;
668                 }
669         }
670
671         /* sort ring by arc length
672          * using a rather bogus insertion sort
673          * butrings will never get too big to matter
674          * */
675         for (i = 0; i < total; i++)
676         {
677                 int j;
678
679                 for (j = i - 1; j >= 0; j--)
680                 {
681                         BArc *arc1, *arc2;
682                         
683                         arc1 = ring[j].arc;
684                         arc2 = ring[j + 1].arc;
685                         
686                         if (arc1->length > arc2->length)
687                         {
688                                 /* swap with smaller */
689                                 RadialArc tmp;
690                                 
691                                 tmp = ring[j + 1];
692                                 ring[j + 1] = ring[j];
693                                 ring[j] = tmp;
694                         }
695                         else
696                         {
697                                 break;
698                         }
699                 }
700         }
701
702         /* Dispatch to specific symmetry tests */
703         first = 0;
704         group = 0;
705         
706         for (i = 1; i < total; i++)
707         {
708                 int dispatch = 0;
709                 int last = i - 1;
710                 
711                 if (fabs(ring[first].arc->length - ring[i].arc->length) > limit)
712                 {
713                         dispatch = 1;
714                 }
715
716                 /* if not dispatching already and on last arc
717                  * Dispatch using current arc as last
718                  * */           
719                 if (dispatch == 0 && i == total - 1)
720                 {
721                         last = i;
722                         dispatch = 1;
723                 } 
724                 
725                 if (dispatch)
726                 {
727                         int sub_total = last - first + 1; 
728
729                         group += 1;
730
731                         if (sub_total == 1)
732                         {
733                                 group -= 1; /* not really a group so decrement */
734                                 /* NOTHING TO DO */
735                         }
736                         else if (sub_total == 2)
737                         {
738                                 BArc *arc1, *arc2;
739                                 BNode *node1, *node2;
740                                 
741                                 arc1 = ring[first].arc;
742                                 arc2 = ring[last].arc;
743                                 
744                                 node1 = BLI_otherNode(arc1, root_node);
745                                 node2 = BLI_otherNode(arc2, root_node);
746                                 
747                                 testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, group);
748                         }
749                         else if (sub_total != total) /* allocate a new sub ring if needed */
750                         {
751                                 RadialArc *sub_ring = MEM_callocN(sizeof(RadialArc) * sub_total, "radial symmetry ring");
752                                 int sub_i;
753                                 
754                                 /* fill in the sub ring */
755                                 for (sub_i = 0; sub_i < sub_total; sub_i++)
756                                 {
757                                         sub_ring[sub_i] = ring[first + sub_i];
758                                 }
759                                 
760                                 testRadialSymmetry(graph, root_node, sub_ring, sub_total, axis, limit, group);
761                         
762                                 MEM_freeN(sub_ring);
763                         }
764                         else if (sub_total == total)
765                         {
766                                 testRadialSymmetry(graph, root_node, ring, total, axis, limit, group);
767                         }
768                         
769                         first = i;
770                 }
771         }
772
773
774         MEM_freeN(ring);
775 }
776
777 static void flagAxialSymmetry(BNode *root_node, BNode *end_node, BArc *arc, int group)
778 {
779         float vec[3];
780         
781         arc->symmetry_group = group;
782         
783         sub_v3_v3v3(vec, end_node->p, root_node->p);
784         
785         if (dot_v3v3(vec, root_node->symmetry_axis) < 0)
786         {
787                 arc->symmetry_flag |= SYM_SIDE_NEGATIVE;
788         }
789         else
790         {
791                 arc->symmetry_flag |= SYM_SIDE_POSITIVE;
792         }
793 }
794
795 static void testAxialSymmetry(BGraph *graph, BNode* root_node, BNode* node1, BNode* node2, BArc* arc1, BArc* arc2, float axis[3], float limit, int group)
796 {
797         float nor[3], vec[3], p[3];
798
799         sub_v3_v3v3(p, node1->p, root_node->p);
800         cross_v3_v3v3(nor, p, axis);
801
802         sub_v3_v3v3(p, root_node->p, node2->p);
803         cross_v3_v3v3(vec, p, axis);
804         add_v3_v3v3(vec, vec, nor);
805         
806         cross_v3_v3v3(nor, vec, axis);
807         
808         if (abs(nor[0]) > abs(nor[1]) && abs(nor[0]) > abs(nor[2]) && nor[0] < 0)
809         {
810                 negate_v3(nor);
811         }
812         else if (abs(nor[1]) > abs(nor[0]) && abs(nor[1]) > abs(nor[2]) && nor[1] < 0)
813         {
814                 negate_v3(nor);
815         }
816         else if (abs(nor[2]) > abs(nor[1]) && abs(nor[2]) > abs(nor[0]) && nor[2] < 0)
817         {
818                 negate_v3(nor);
819         }
820         
821         /* mirror node2 along axis */
822         VECCOPY(p, node2->p);
823         BLI_mirrorAlongAxis(p, root_node->p, nor);
824         
825         /* check if it's within limit before continuing */
826         if (len_v3v3(node1->p, p) <= limit)
827         {
828                 /* mark node as symmetric physically */
829                 VECCOPY(root_node->symmetry_axis, nor);
830                 root_node->symmetry_flag |= SYM_PHYSICAL;
831                 root_node->symmetry_flag |= SYM_AXIAL;
832
833                 /* flag side on arcs */
834                 flagAxialSymmetry(root_node, node1, arc1, group);
835                 flagAxialSymmetry(root_node, node2, arc2, group);
836                 
837                 if (graph->axial_symmetry)
838                 {
839                         graph->axial_symmetry(root_node, node1, node2, arc1, arc2);
840                 }
841         }
842         else
843         {
844                 /* NOT SYMMETRIC */
845         }
846 }
847
848 static void handleAxialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit)
849 {
850         BArc *arc1 = NULL, *arc2 = NULL;
851         BNode *node1 = NULL, *node2 = NULL;
852         int i;
853         
854         /* mark topological symmetry */
855         root_node->symmetry_flag |= SYM_TOPOLOGICAL;
856
857         for (i = 0; i < root_node->degree; i++)
858         {
859                 BArc *connectedArc = root_node->arcs[i];
860                 
861                 /* depth is store as a negative in flag. symmetry level is positive */
862                 if (connectedArc->symmetry_level == -depth)
863                 {
864                         if (arc1 == NULL)
865                         {
866                                 arc1 = connectedArc;
867                                 node1 = BLI_otherNode(arc1, root_node);
868                         }
869                         else
870                         {
871                                 arc2 = connectedArc;
872                                 node2 = BLI_otherNode(arc2, root_node);
873                                 break; /* Can stop now, the two arcs have been found */
874                         }
875                 }
876         }
877         
878         /* shouldn't happen, but just to be sure */
879         if (node1 == NULL || node2 == NULL)
880         {
881                 return;
882         }
883         
884         testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, 1);
885 }
886
887 static void markdownSecondarySymmetry(BGraph *graph, BNode *node, int depth, int level, float limit)
888 {
889         float axis[3] = {0, 0, 0};
890         int count = 0;
891         int i;
892         
893         /* count the number of branches in this symmetry group
894          * and determinte the axis of symmetry
895          *  */  
896         for (i = 0; i < node->degree; i++)
897         {
898                 BArc *connectedArc = node->arcs[i];
899                 
900                 /* depth is store as a negative in flag. symmetry level is positive */
901                 if (connectedArc->symmetry_level == -depth)
902                 {
903                         count++;
904                 }
905                 /* If arc is on the axis */
906                 else if (connectedArc->symmetry_level == level)
907                 {
908                         add_v3_v3v3(axis, axis, connectedArc->head->p);
909                         sub_v3_v3v3(axis, axis, connectedArc->tail->p);
910                 }
911         }
912
913         normalize_v3(axis);
914
915         /* Split between axial and radial symmetry */
916         if (count == 2)
917         {
918                 handleAxialSymmetry(graph, node, depth, axis, limit);
919         }
920         else
921         {
922                 handleRadialSymmetry(graph, node, depth, axis, limit);
923         }
924                 
925         /* markdown secondary symetries */      
926         for (i = 0; i < node->degree; i++)
927         {
928                 BArc *connectedArc = node->arcs[i];
929                 
930                 if (connectedArc->symmetry_level == -depth)
931                 {
932                         /* markdown symmetry for branches corresponding to the depth */
933                         markdownSymmetryArc(graph, connectedArc, node, level + 1, limit);
934                 }
935         }
936 }
937
938 static void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit)
939 {
940         int i;
941
942         /* if arc is null, we start straight from a node */     
943         if (arc)
944         {
945                 arc->symmetry_level = level;
946                 
947                 node = BLI_otherNode(arc, node);
948         }
949         
950         for (i = 0; i < node->degree; i++)
951         {
952                 BArc *connectedArc = node->arcs[i];
953                 
954                 if (connectedArc != arc)
955                 {
956                         BNode *connectedNode = BLI_otherNode(connectedArc, node);
957                         
958                         /* symmetry level is positive value, negative values is subtree depth */
959                         connectedArc->symmetry_level = -BLI_subtreeShape(graph, connectedNode, connectedArc, 0);
960                 }
961         }
962
963         arc = NULL;
964
965         for (i = 0; i < node->degree; i++)
966         {
967                 int issymmetryAxis = 0;
968                 BArc *connectedArc = node->arcs[i];
969                 
970                 /* only arcs not already marked as symetric */
971                 if (connectedArc->symmetry_level < 0)
972                 {
973                         int j;
974                         
975                         /* true by default */
976                         issymmetryAxis = 1;
977                         
978                         for (j = 0; j < node->degree; j++)
979                         {
980                                 BArc *otherArc = node->arcs[j];
981                                 
982                                 /* different arc, same depth */
983                                 if (otherArc != connectedArc && otherArc->symmetry_level == connectedArc->symmetry_level)
984                                 {
985                                         /* not on the symmetry axis */
986                                         issymmetryAxis = 0;
987                                         break;
988                                 } 
989                         }
990                 }
991                 
992                 /* arc could be on the symmetry axis */
993                 if (issymmetryAxis == 1)
994                 {
995                         /* no arc as been marked previously, keep this one */
996                         if (arc == NULL)
997                         {
998                                 arc = connectedArc;
999                         }
1000                         else if (connectedArc->symmetry_level < arc->symmetry_level)
1001                         {
1002                                 /* go with more complex subtree as symmetry arc */
1003                                 arc = connectedArc;
1004                         }
1005                 }
1006         }
1007         
1008         /* go down the arc continuing the symmetry axis */
1009         if (arc)
1010         {
1011                 markdownSymmetryArc(graph, arc, node, level, limit);
1012         }
1013
1014         
1015         /* secondary symmetry */
1016         for (i = 0; i < node->degree; i++)
1017         {
1018                 BArc *connectedArc = node->arcs[i];
1019                 
1020                 /* only arcs not already marked as symetric and is not the next arc on the symmetry axis */
1021                 if (connectedArc->symmetry_level < 0)
1022                 {
1023                         /* subtree depth is store as a negative value in the symmetry */
1024                         markdownSecondarySymmetry(graph, node, -connectedArc->symmetry_level, level, limit);
1025                 }
1026         }
1027 }
1028
1029 void BLI_markdownSymmetry(BGraph *graph, BNode *root_node, float limit)
1030 {
1031         BNode *node;
1032         BArc *arc;
1033         
1034         if (root_node == NULL)
1035         {
1036                 return;
1037         }
1038         
1039         if (BLI_isGraphCyclic(graph))
1040         {
1041                 return;
1042         }
1043         
1044         /* mark down all arcs as non-symetric */
1045         BLI_flagArcs(graph, 0);
1046         
1047         /* mark down all nodes as not on the symmetry axis */
1048         BLI_flagNodes(graph, 0);
1049
1050         node = root_node;
1051         
1052         /* sanity check REMOVE ME */
1053         if (node->degree > 0)
1054         {
1055                 arc = node->arcs[0];
1056                 
1057                 if (node->degree == 1)
1058                 {
1059                         markdownSymmetryArc(graph, arc, node, 1, limit);
1060                 }
1061                 else
1062                 {
1063                         markdownSymmetryArc(graph, NULL, node, 1, limit);
1064                 }
1065                 
1066
1067
1068                 /* mark down non-symetric arcs */
1069                 for (arc = graph->arcs.first; arc; arc = arc->next)
1070                 {
1071                         if (arc->symmetry_level < 0)
1072                         {
1073                                 arc->symmetry_level = 0;
1074                         }
1075                         else
1076                         {
1077                                 /* mark down nodes with the lowest level symmetry axis */
1078                                 if (arc->head->symmetry_level == 0 || arc->head->symmetry_level > arc->symmetry_level)
1079                                 {
1080                                         arc->head->symmetry_level = arc->symmetry_level;
1081                                 }
1082                                 if (arc->tail->symmetry_level == 0 || arc->tail->symmetry_level > arc->symmetry_level)
1083                                 {
1084                                         arc->tail->symmetry_level = arc->symmetry_level;
1085                                 }
1086                         }
1087                 }
1088         }
1089 }
1090
1091 void* IT_head(void* arg)
1092 {
1093         BArcIterator *iter = (BArcIterator*)arg; 
1094         return iter->head(iter);
1095 }
1096
1097 void* IT_tail(void* arg)
1098 {
1099         BArcIterator *iter = (BArcIterator*)arg; 
1100         return iter->tail(iter); 
1101 }
1102
1103 void* IT_peek(void* arg, int n)
1104 {
1105         BArcIterator *iter = (BArcIterator*)arg;
1106         
1107         if (iter->index + n < 0)
1108         {
1109                 return iter->head(iter);
1110         }
1111         else if (iter->index + n >= iter->length)
1112         {
1113                 return iter->tail(iter);
1114         }
1115         else
1116         {
1117                 return iter->peek(iter, n);
1118         }
1119 }
1120
1121 void* IT_next(void* arg)
1122 {
1123         BArcIterator *iter = (BArcIterator*)arg; 
1124         return iter->next(iter);
1125 }
1126
1127 void* IT_nextN(void* arg, int n)
1128 {
1129         BArcIterator *iter = (BArcIterator*)arg; 
1130         return iter->nextN(iter, n);
1131 }
1132
1133 void* IT_previous(void* arg)
1134 {
1135         BArcIterator *iter = (BArcIterator*)arg; 
1136         return iter->previous(iter);
1137 }
1138
1139 int   IT_stopped(void* arg)
1140 {
1141         BArcIterator *iter = (BArcIterator*)arg; 
1142         return iter->stopped(iter);
1143 }