Patch from Banlu Kemiyatorn
[blender.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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, 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_arithb.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 && VecLenf(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;
283         
284         for(node = graph->nodes.first; node; node = node->next)
285         {
286                 float distance = VecLenf(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 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         BNode *test_node;
469         
470         BLI_flagNodes(graph, 0);
471         return subtreeShape(node, rootArc, include_root);
472 }
473
474 float BLI_subtreeLength(BNode *node)
475 {
476         float length = 0;
477         int i;
478
479         node->flag = 0; /* flag node as visited */
480
481         for(i = 0; i < node->degree; i++)
482         {
483                 BArc *arc = node->arcs[i];
484                 BNode *other_node = BLI_otherNode(arc, node);
485                 
486                 if (other_node->flag != 0)
487                 {
488                         float subgraph_length = arc->length + BLI_subtreeLength(other_node); 
489                         length = MAX2(length, subgraph_length);
490                 }
491         }
492         
493         return length;
494 }
495
496 void BLI_calcGraphLength(BGraph *graph)
497 {
498         float length = 0;
499         int nb_subgraphs;
500         int i;
501         
502         nb_subgraphs = BLI_FlagSubgraphs(graph);
503         
504         for (i = 1; i <= nb_subgraphs; i++)
505         {
506                 BNode *node;
507                 
508                 for (node = graph->nodes.first; node; node = node->next)
509                 {
510                         /* start on an external node  of the subgraph */
511                         if (node->subgraph_index == i && node->degree == 1)
512                         {
513                                 float subgraph_length = BLI_subtreeLength(node);
514                                 length = MAX2(length, subgraph_length);
515                                 break;
516                         }
517                 }
518         }
519         
520         graph->length = length;
521 }
522
523 /********************************* SYMMETRY DETECTION **************************************************/
524
525 void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit);
526
527 void BLI_mirrorAlongAxis(float v[3], float center[3], float axis[3])
528 {
529         float dv[3], pv[3];
530         
531         VecSubf(dv, v, center);
532         Projf(pv, dv, axis);
533         VecMulf(pv, -2);
534         VecAddf(v, v, pv);
535 }
536
537 static void testRadialSymmetry(BGraph *graph, BNode* root_node, RadialArc* ring, int total, float axis[3], float limit, int group)
538 {
539         int symmetric = 1;
540         int i;
541         
542         /* sort ring by angle */
543         for (i = 0; i < total - 1; i++)
544         {
545                 float minAngle = FLT_MAX;
546                 int minIndex = -1;
547                 int j;
548
549                 for (j = i + 1; j < total; j++)
550                 {
551                         float angle = Inpf(ring[i].n, ring[j].n);
552
553                         /* map negative values to 1..2 */
554                         if (angle < 0)
555                         {
556                                 angle = 1 - angle;
557                         }
558
559                         if (angle < minAngle)
560                         {
561                                 minIndex = j;
562                                 minAngle = angle;
563                         }
564                 }
565
566                 /* swap if needed */
567                 if (minIndex != i + 1)
568                 {
569                         RadialArc tmp;
570                         tmp = ring[i + 1];
571                         ring[i + 1] = ring[minIndex];
572                         ring[minIndex] = tmp;
573                 }
574         }
575
576         for (i = 0; i < total && symmetric; i++)
577         {
578                 BNode *node1, *node2;
579                 float tangent[3];
580                 float normal[3];
581                 float p[3];
582                 int j = (i + 1) % total; /* next arc in the circular list */
583
584                 VecAddf(tangent, ring[i].n, ring[j].n);
585                 Crossf(normal, tangent, axis);
586                 
587                 node1 = BLI_otherNode(ring[i].arc, root_node);
588                 node2 = BLI_otherNode(ring[j].arc, root_node);
589
590                 VECCOPY(p, node2->p);
591                 BLI_mirrorAlongAxis(p, root_node->p, normal);
592                 
593                 /* check if it's within limit before continuing */
594                 if (VecLenf(node1->p, p) > limit)
595                 {
596                         symmetric = 0;
597                 }
598
599         }
600
601         if (symmetric)
602         {
603                 /* mark node as symmetric physically */
604                 VECCOPY(root_node->symmetry_axis, axis);
605                 root_node->symmetry_flag |= SYM_PHYSICAL;
606                 root_node->symmetry_flag |= SYM_RADIAL;
607                 
608                 /* FLAG SYMMETRY GROUP */
609                 for (i = 0; i < total; i++)
610                 {
611                         ring[i].arc->symmetry_group = group;
612                         ring[i].arc->symmetry_flag = SYM_SIDE_RADIAL + i;
613                 }
614
615                 if (graph->radial_symmetry)
616                 {
617                         graph->radial_symmetry(root_node, ring, total);
618                 }
619         }
620 }
621
622 static void handleRadialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit)
623 {
624         RadialArc *ring = NULL;
625         RadialArc *unit;
626         int total = 0;
627         int group;
628         int first;
629         int i;
630
631         /* mark topological symmetry */
632         root_node->symmetry_flag |= SYM_TOPOLOGICAL;
633
634         /* total the number of arcs in the symmetry ring */
635         for (i = 0; i < root_node->degree; i++)
636         {
637                 BArc *connectedArc = root_node->arcs[i];
638                 
639                 /* depth is store as a negative in flag. symmetry level is positive */
640                 if (connectedArc->symmetry_level == -depth)
641                 {
642                         total++;
643                 }
644         }
645
646         ring = MEM_callocN(sizeof(RadialArc) * total, "radial symmetry ring");
647         unit = ring;
648
649         /* fill in the ring */
650         for (unit = ring, i = 0; i < root_node->degree; i++)
651         {
652                 BArc *connectedArc = root_node->arcs[i];
653                 
654                 /* depth is store as a negative in flag. symmetry level is positive */
655                 if (connectedArc->symmetry_level == -depth)
656                 {
657                         BNode *otherNode = BLI_otherNode(connectedArc, root_node);
658                         float vec[3];
659
660                         unit->arc = connectedArc;
661
662                         /* project the node to node vector on the symmetry plane */
663                         VecSubf(unit->n, otherNode->p, root_node->p);
664                         Projf(vec, unit->n, axis);
665                         VecSubf(unit->n, unit->n, vec);
666
667                         Normalize(unit->n);
668
669                         unit++;
670                 }
671         }
672
673         /* sort ring by arc length
674          * using a rather bogus insertion sort
675          * butrings will never get too big to matter
676          * */
677         for (i = 0; i < total; i++)
678         {
679                 int j;
680
681                 for (j = i - 1; j >= 0; j--)
682                 {
683                         BArc *arc1, *arc2;
684                         
685                         arc1 = ring[j].arc;
686                         arc2 = ring[j + 1].arc;
687                         
688                         if (arc1->length > arc2->length)
689                         {
690                                 /* swap with smaller */
691                                 RadialArc tmp;
692                                 
693                                 tmp = ring[j + 1];
694                                 ring[j + 1] = ring[j];
695                                 ring[j] = tmp;
696                         }
697                         else
698                         {
699                                 break;
700                         }
701                 }
702         }
703
704         /* Dispatch to specific symmetry tests */
705         first = 0;
706         group = 0;
707         
708         for (i = 1; i < total; i++)
709         {
710                 int dispatch = 0;
711                 int last = i - 1;
712                 
713                 if (fabs(ring[first].arc->length - ring[i].arc->length) > limit)
714                 {
715                         dispatch = 1;
716                 }
717
718                 /* if not dispatching already and on last arc
719                  * Dispatch using current arc as last
720                  * */           
721                 if (dispatch == 0 && i == total - 1)
722                 {
723                         last = i;
724                         dispatch = 1;
725                 } 
726                 
727                 if (dispatch)
728                 {
729                         int sub_total = last - first + 1; 
730
731                         group += 1;
732
733                         if (sub_total == 1)
734                         {
735                                 group -= 1; /* not really a group so decrement */
736                                 /* NOTHING TO DO */
737                         }
738                         else if (sub_total == 2)
739                         {
740                                 BArc *arc1, *arc2;
741                                 BNode *node1, *node2;
742                                 
743                                 arc1 = ring[first].arc;
744                                 arc2 = ring[last].arc;
745                                 
746                                 node1 = BLI_otherNode(arc1, root_node);
747                                 node2 = BLI_otherNode(arc2, root_node);
748                                 
749                                 testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, group);
750                         }
751                         else if (sub_total != total) /* allocate a new sub ring if needed */
752                         {
753                                 RadialArc *sub_ring = MEM_callocN(sizeof(RadialArc) * sub_total, "radial symmetry ring");
754                                 int sub_i;
755                                 
756                                 /* fill in the sub ring */
757                                 for (sub_i = 0; sub_i < sub_total; sub_i++)
758                                 {
759                                         sub_ring[sub_i] = ring[first + sub_i];
760                                 }
761                                 
762                                 testRadialSymmetry(graph, root_node, sub_ring, sub_total, axis, limit, group);
763                         
764                                 MEM_freeN(sub_ring);
765                         }
766                         else if (sub_total == total)
767                         {
768                                 testRadialSymmetry(graph, root_node, ring, total, axis, limit, group);
769                         }
770                         
771                         first = i;
772                 }
773         }
774
775
776         MEM_freeN(ring);
777 }
778
779 static void flagAxialSymmetry(BNode *root_node, BNode *end_node, BArc *arc, int group)
780 {
781         float vec[3];
782         
783         arc->symmetry_group = group;
784         
785         VecSubf(vec, end_node->p, root_node->p);
786         
787         if (Inpf(vec, root_node->symmetry_axis) < 0)
788         {
789                 arc->symmetry_flag |= SYM_SIDE_NEGATIVE;
790         }
791         else
792         {
793                 arc->symmetry_flag |= SYM_SIDE_POSITIVE;
794         }
795 }
796
797 static void testAxialSymmetry(BGraph *graph, BNode* root_node, BNode* node1, BNode* node2, BArc* arc1, BArc* arc2, float axis[3], float limit, int group)
798 {
799         float nor[3], vec[3], p[3];
800
801         VecSubf(p, node1->p, root_node->p);
802         Crossf(nor, p, axis);
803
804         VecSubf(p, root_node->p, node2->p);
805         Crossf(vec, p, axis);
806         VecAddf(vec, vec, nor);
807         
808         Crossf(nor, vec, axis);
809         
810         if (abs(nor[0]) > abs(nor[1]) && abs(nor[0]) > abs(nor[2]) && nor[0] < 0)
811         {
812                 VecNegf(nor);
813         }
814         else if (abs(nor[1]) > abs(nor[0]) && abs(nor[1]) > abs(nor[2]) && nor[1] < 0)
815         {
816                 VecNegf(nor);
817         }
818         else if (abs(nor[2]) > abs(nor[1]) && abs(nor[2]) > abs(nor[0]) && nor[2] < 0)
819         {
820                 VecNegf(nor);
821         }
822         
823         /* mirror node2 along axis */
824         VECCOPY(p, node2->p);
825         BLI_mirrorAlongAxis(p, root_node->p, nor);
826         
827         /* check if it's within limit before continuing */
828         if (VecLenf(node1->p, p) <= limit)
829         {
830                 /* mark node as symmetric physically */
831                 VECCOPY(root_node->symmetry_axis, nor);
832                 root_node->symmetry_flag |= SYM_PHYSICAL;
833                 root_node->symmetry_flag |= SYM_AXIAL;
834
835                 /* flag side on arcs */
836                 flagAxialSymmetry(root_node, node1, arc1, group);
837                 flagAxialSymmetry(root_node, node2, arc2, group);
838                 
839                 if (graph->axial_symmetry)
840                 {
841                         graph->axial_symmetry(root_node, node1, node2, arc1, arc2);
842                 }
843         }
844         else
845         {
846                 /* NOT SYMMETRIC */
847         }
848 }
849
850 static void handleAxialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit)
851 {
852         BArc *arc1 = NULL, *arc2 = NULL;
853         BNode *node1 = NULL, *node2 = NULL;
854         int i;
855         
856         /* mark topological symmetry */
857         root_node->symmetry_flag |= SYM_TOPOLOGICAL;
858
859         for (i = 0; i < root_node->degree; i++)
860         {
861                 BArc *connectedArc = root_node->arcs[i];
862                 
863                 /* depth is store as a negative in flag. symmetry level is positive */
864                 if (connectedArc->symmetry_level == -depth)
865                 {
866                         if (arc1 == NULL)
867                         {
868                                 arc1 = connectedArc;
869                                 node1 = BLI_otherNode(arc1, root_node);
870                         }
871                         else
872                         {
873                                 arc2 = connectedArc;
874                                 node2 = BLI_otherNode(arc2, root_node);
875                                 break; /* Can stop now, the two arcs have been found */
876                         }
877                 }
878         }
879         
880         /* shouldn't happen, but just to be sure */
881         if (node1 == NULL || node2 == NULL)
882         {
883                 return;
884         }
885         
886         testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, 1);
887 }
888
889 static void markdownSecondarySymmetry(BGraph *graph, BNode *node, int depth, int level, float limit)
890 {
891         float axis[3] = {0, 0, 0};
892         int count = 0;
893         int i;
894         
895         /* count the number of branches in this symmetry group
896          * and determinte the axis of symmetry
897          *  */  
898         for (i = 0; i < node->degree; i++)
899         {
900                 BArc *connectedArc = node->arcs[i];
901                 
902                 /* depth is store as a negative in flag. symmetry level is positive */
903                 if (connectedArc->symmetry_level == -depth)
904                 {
905                         count++;
906                 }
907                 /* If arc is on the axis */
908                 else if (connectedArc->symmetry_level == level)
909                 {
910                         VecAddf(axis, axis, connectedArc->head->p);
911                         VecSubf(axis, axis, connectedArc->tail->p);
912                 }
913         }
914
915         Normalize(axis);
916
917         /* Split between axial and radial symmetry */
918         if (count == 2)
919         {
920                 handleAxialSymmetry(graph, node, depth, axis, limit);
921         }
922         else
923         {
924                 handleRadialSymmetry(graph, node, depth, axis, limit);
925         }
926                 
927         /* markdown secondary symetries */      
928         for (i = 0; i < node->degree; i++)
929         {
930                 BArc *connectedArc = node->arcs[i];
931                 
932                 if (connectedArc->symmetry_level == -depth)
933                 {
934                         /* markdown symmetry for branches corresponding to the depth */
935                         markdownSymmetryArc(graph, connectedArc, node, level + 1, limit);
936                 }
937         }
938 }
939
940 void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit)
941 {
942         int i;
943
944         /* if arc is null, we start straight from a node */     
945         if (arc)
946         {
947                 arc->symmetry_level = level;
948                 
949                 node = BLI_otherNode(arc, node);
950         }
951         
952         for (i = 0; i < node->degree; i++)
953         {
954                 BArc *connectedArc = node->arcs[i];
955                 
956                 if (connectedArc != arc)
957                 {
958                         BNode *connectedNode = BLI_otherNode(connectedArc, node);
959                         
960                         /* symmetry level is positive value, negative values is subtree depth */
961                         connectedArc->symmetry_level = -BLI_subtreeShape(graph, connectedNode, connectedArc, 0);
962                 }
963         }
964
965         arc = NULL;
966
967         for (i = 0; i < node->degree; i++)
968         {
969                 int issymmetryAxis = 0;
970                 BArc *connectedArc = node->arcs[i];
971                 
972                 /* only arcs not already marked as symetric */
973                 if (connectedArc->symmetry_level < 0)
974                 {
975                         int j;
976                         
977                         /* true by default */
978                         issymmetryAxis = 1;
979                         
980                         for (j = 0; j < node->degree; j++)
981                         {
982                                 BArc *otherArc = node->arcs[j];
983                                 
984                                 /* different arc, same depth */
985                                 if (otherArc != connectedArc && otherArc->symmetry_level == connectedArc->symmetry_level)
986                                 {
987                                         /* not on the symmetry axis */
988                                         issymmetryAxis = 0;
989                                         break;
990                                 } 
991                         }
992                 }
993                 
994                 /* arc could be on the symmetry axis */
995                 if (issymmetryAxis == 1)
996                 {
997                         /* no arc as been marked previously, keep this one */
998                         if (arc == NULL)
999                         {
1000                                 arc = connectedArc;
1001                         }
1002                         else if (connectedArc->symmetry_level < arc->symmetry_level)
1003                         {
1004                                 /* go with more complex subtree as symmetry arc */
1005                                 arc = connectedArc;
1006                         }
1007                 }
1008         }
1009         
1010         /* go down the arc continuing the symmetry axis */
1011         if (arc)
1012         {
1013                 markdownSymmetryArc(graph, arc, node, level, limit);
1014         }
1015
1016         
1017         /* secondary symmetry */
1018         for (i = 0; i < node->degree; i++)
1019         {
1020                 BArc *connectedArc = node->arcs[i];
1021                 
1022                 /* only arcs not already marked as symetric and is not the next arc on the symmetry axis */
1023                 if (connectedArc->symmetry_level < 0)
1024                 {
1025                         /* subtree depth is store as a negative value in the symmetry */
1026                         markdownSecondarySymmetry(graph, node, -connectedArc->symmetry_level, level, limit);
1027                 }
1028         }
1029 }
1030
1031 void BLI_markdownSymmetry(BGraph *graph, BNode *root_node, float limit)
1032 {
1033         BNode *node;
1034         BArc *arc;
1035         
1036         if (BLI_isGraphCyclic(graph))
1037         {
1038                 return;
1039         }
1040         
1041         /* mark down all arcs as non-symetric */
1042         BLI_flagArcs(graph, 0);
1043         
1044         /* mark down all nodes as not on the symmetry axis */
1045         BLI_flagNodes(graph, 0);
1046
1047         node = root_node;
1048         
1049         /* sanity check REMOVE ME */
1050         if (node->degree > 0)
1051         {
1052                 arc = node->arcs[0];
1053                 
1054                 if (node->degree == 1)
1055                 {
1056                         markdownSymmetryArc(graph, arc, node, 1, limit);
1057                 }
1058                 else
1059                 {
1060                         markdownSymmetryArc(graph, NULL, node, 1, limit);
1061                 }
1062                 
1063
1064
1065                 /* mark down non-symetric arcs */
1066                 for (arc = graph->arcs.first; arc; arc = arc->next)
1067                 {
1068                         if (arc->symmetry_level < 0)
1069                         {
1070                                 arc->symmetry_level = 0;
1071                         }
1072                         else
1073                         {
1074                                 /* mark down nodes with the lowest level symmetry axis */
1075                                 if (arc->head->symmetry_level == 0 || arc->head->symmetry_level > arc->symmetry_level)
1076                                 {
1077                                         arc->head->symmetry_level = arc->symmetry_level;
1078                                 }
1079                                 if (arc->tail->symmetry_level == 0 || arc->tail->symmetry_level > arc->symmetry_level)
1080                                 {
1081                                         arc->tail->symmetry_level = arc->symmetry_level;
1082                                 }
1083                         }
1084                 }
1085         }
1086 }
1087