Cleaned up radiosity mesh adding code, was using ugly pointer tricks.
[blender.git] / source / blender / radiosity / intern / source / radnode.c
1 /* ***************************************
2  *
3  * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version. The Blender
9  * Foundation also sells licenses for use in proprietary software under
10  * the Blender License.  See http://www.blender.org/BL/ for information
11  * about this.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software Foundation,
20  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
21  *
22  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
23  * All rights reserved.
24  *
25  * The Original Code is: all of this file.
26  *
27  * Contributor(s): none yet.
28  *
29  * ***** END GPL/BL DUAL LICENSE BLOCK *****
30
31
32
33     node.c      nov/dec 1992
34                         may 1999
35
36     $Id$
37
38  *************************************** */
39
40 #include <stdio.h>
41 #include <math.h>
42 #include <string.h>
43
44 #include "MEM_guardedalloc.h"
45
46 #include "BLI_blenlib.h"
47 #include "BLI_arithb.h"
48
49 #include "BKE_global.h"
50 #include "BKE_main.h"
51
52 #include "BIF_toolbox.h"
53
54 #include "radio.h"
55
56 #ifdef HAVE_CONFIG_H
57 #include <config.h>
58 #endif
59
60 /* locals */
61 void *malloc_fast(int size);
62 void *calloc_fast(int size);
63 void free_fast(void *poin, int siz);
64 void deleteTriNodes(RNode *node);
65 /* lower because of local type define */
66 /*  void check_mallocgroup(MallocGroup *mg); */
67
68
69 /* ********** fastmalloc ************** */
70
71 #define MAL_GROUPSIZE   256
72 #define MAL_AVAILABLE   1
73 #define MAL_FULL                2
74
75
76
77
78 ListBase MallocBase= {0, 0};
79 int totfastmem= 0;
80
81 typedef struct MallocGroup {
82         struct MallocGroup *next, *prev;
83         short size, flag;
84         short curfree, tot;
85         char flags[MAL_GROUPSIZE];
86         char *data;
87 } MallocGroup;
88
89 /* one more local */
90 void check_mallocgroup(MallocGroup *mg);
91
92 void check_mallocgroup(MallocGroup *mg)
93 {
94         int a;
95         char *cp;
96         
97         if(mg->tot==MAL_GROUPSIZE) {
98                 mg->flag= MAL_FULL;
99                 return;
100         }
101         
102         cp= mg->flags;
103         
104         if(mg->curfree<MAL_GROUPSIZE-1) {
105                 if(cp[mg->curfree+1]==0) {
106                         mg->curfree++;
107                         return;
108                 }
109         }
110         if(mg->curfree>0) {
111                 if(cp[mg->curfree-1]==0) {
112                         mg->curfree--;
113                         return;
114                 }
115         }
116
117         for(a=0; a<MAL_GROUPSIZE; a++) {
118                 if(cp[a]==0) {
119                         mg->curfree= a;
120                         return;
121                 }
122         }
123         printf("fastmalloc: shouldnt be here\n");
124 }
125
126 void *malloc_fast(int size)
127 {
128         MallocGroup *mg;
129         void *retval;
130         
131         mg= MallocBase.last;
132         while(mg) {
133                 if(mg->size==size) {
134                         if(mg->flag & MAL_AVAILABLE) {
135                                 mg->flags[mg->curfree]= 1;
136                                 mg->tot++;
137                                 retval= mg->data+mg->curfree*mg->size;
138                                 check_mallocgroup(mg);
139                                 return retval;
140                         }
141                 }
142                 mg= mg->prev;
143         }
144
145         /* no free block found */
146         mg= MEM_callocN(sizeof(MallocGroup), "mallocgroup");
147         BLI_addtail(&MallocBase, mg);
148         mg->data= MEM_mallocN(MAL_GROUPSIZE*size, "mallocgroupdata");
149         mg->flag= MAL_AVAILABLE;
150         mg->flags[0]= 1;
151         mg->curfree= 1;
152         mg->size= size;
153         mg->tot= 1;
154         
155         totfastmem+= sizeof(MallocGroup)+MAL_GROUPSIZE*size;
156
157         return mg->data;
158 }
159
160 void *calloc_fast(int size)
161 {
162         void *poin;
163         
164         poin= malloc_fast(size);
165         memset(poin, 0, size);
166         
167         return poin;
168 }
169
170 void free_fast(void *poin, int size)
171 {
172         MallocGroup *mg;
173         long val;
174
175         mg= MallocBase.last;
176         while(mg) {
177                 if(mg->size==size) {
178                         if( ((long)poin) >= ((long)mg->data) ) {
179                                 if( ((long)poin) < ((long)(mg->data+MAL_GROUPSIZE*size)) ) {
180                                         val= ((long)poin) - ((long)mg->data);
181                                         val/= size;
182                                         mg->curfree= val;
183                                         mg->flags[val]= 0;
184                                         mg->flag= MAL_AVAILABLE;
185                                         
186                                         mg->tot--;
187                                         if(mg->tot==0) {
188                                                 BLI_remlink(&MallocBase, mg);
189                                                 MEM_freeN(mg->data);
190                                                 MEM_freeN(mg);
191                                                 totfastmem-= sizeof(MallocGroup)+MAL_GROUPSIZE*size;
192                                         }
193                                         return;
194                                 }
195                         }
196                 }
197                 mg= mg->prev;
198         }
199         printf("fast free: pointer not in memlist %p size %d\n",
200                    poin, size);
201 }
202
203 /* security: only one function in a time can use it */
204 static char *fastmallocstr= 0;
205
206 void free_fastAll()
207 {
208         MallocGroup *mg;
209         
210         mg= MallocBase.first;
211         while(mg) {
212                 BLI_remlink(&MallocBase, mg);
213                 MEM_freeN(mg->data);
214                 MEM_freeN(mg);
215                 mg= MallocBase.first;
216         }
217         totfastmem= 0;
218         fastmallocstr= 0;
219 }
220
221 void start_fastmalloc(char *str)
222 {
223         if(fastmallocstr) {
224                 error("Fastmalloc in use: %s", fastmallocstr);
225                 return;
226         }
227         fastmallocstr= str;
228 }
229
230 /* **************************************** */
231
232 float nodelimit;
233
234 void setnodelimit(float limit)
235 {
236         nodelimit= limit;
237
238 }
239
240 /* ************  memory management ***********  */
241
242 int Ntotvert=0, Ntotnode=0, Ntotpatch=0;
243
244 float *mallocVert()
245 {
246         Ntotvert++;
247         return (float *)malloc_fast(16);
248 }
249
250 float *callocVert()
251 {
252         Ntotvert++;
253         return (float *)calloc_fast(16);
254 }
255
256 void freeVert(float *vert)
257 {
258         free_fast(vert, 16);
259         Ntotvert--;
260 }
261
262 int totalRadVert()
263 {
264         return Ntotvert;
265 }
266
267 RNode *mallocNode()
268 {
269         Ntotnode++;
270         return (RNode *)malloc_fast(sizeof(RNode));
271 }
272
273 RNode *callocNode()
274 {
275         Ntotnode++;
276         return (RNode *)calloc_fast(sizeof(RNode));
277 }
278
279 void freeNode(RNode *node)
280 {
281         free_fast(node, sizeof(RNode));
282         Ntotnode--;
283 }
284
285 void freeNode_recurs(RNode *node)
286 {
287
288         if(node->down1) {
289                 freeNode_recurs(node->down1);
290                 freeNode_recurs(node->down2);
291         }
292
293         node->down1= node->down2= 0;
294         freeNode(node);
295
296 }
297
298 RPatch *mallocPatch()
299 {
300         Ntotpatch++;
301         return (RPatch *)malloc_fast(sizeof(RPatch));
302 }
303
304 RPatch *callocPatch()
305 {
306         Ntotpatch++;
307         return (RPatch *)calloc_fast(sizeof(RPatch));
308 }
309
310 void freePatch(RPatch *patch)
311 {
312         free_fast(patch, sizeof(RPatch));
313         Ntotpatch--;
314 }
315
316 /* ************  SUBDIVIDE  ***********  */
317
318
319 void replaceAllNode(RNode *neighb, RNode *newn)
320 {
321         /* changes from all neighbours the edgepointers that point to newn->up in new */
322         int ok= 0;
323         
324         
325         if(neighb==0) return;
326         if(newn->up==0) return;
327         
328         if(neighb->ed1==newn->up) {
329                 neighb->ed1= newn;
330                 ok= 1;
331         }
332         else if(neighb->ed2==newn->up) {
333                 neighb->ed2= newn;
334                 ok= 1;
335         }
336         else if(neighb->ed3==newn->up) {
337                 neighb->ed3= newn;
338                 ok= 1;
339         }
340         else if(neighb->ed4==newn->up) {
341                 neighb->ed4= newn;
342                 ok= 1;
343         }
344         
345         if(ok && neighb->down1) {
346                 replaceAllNode(neighb->down1, newn);
347                 replaceAllNode(neighb->down2, newn);
348         }
349 }
350
351 void replaceAllNodeInv(RNode *neighb, RNode *old)
352 {
353         /* changes from all neighbours the edgepointers that point to old in old->up */
354         if(neighb==0) return;
355         if(old->up==0) return;
356         
357         if(neighb->ed1==old) {
358                 neighb->ed1= old->up;
359         }
360         else if(neighb->ed2==old) {
361                 neighb->ed2= old->up;
362         }
363         else if(neighb->ed3==old) {
364                 neighb->ed3= old->up;
365         }
366         else if(neighb->ed4==old) {
367                 neighb->ed4= old->up;
368         }
369         
370         if(neighb->down1) {
371                 replaceAllNodeInv(neighb->down1, old);
372                 replaceAllNodeInv(neighb->down2, old);
373         }
374 }
375
376 void replaceAllNodeUp(RNode *neighb, RNode *old)
377 {
378         /* changes from all neighbours the edgepointers that point to old in old->up */
379         if(neighb==0) return;
380         if(old->up==0) return;
381         neighb= neighb->up;
382         if(neighb==0) return;
383         
384         if(neighb->ed1==old) {
385                 neighb->ed1= old->up;
386         }
387         else if(neighb->ed2==old) {
388                 neighb->ed2= old->up;
389         }
390         else if(neighb->ed3==old) {
391                 neighb->ed3= old->up;
392         }
393         else if(neighb->ed4==old) {
394                 neighb->ed4= old->up;
395         }
396         
397         if(neighb->up) {
398                 replaceAllNodeUp(neighb, old);
399         }
400 }
401
402
403 void replaceTestNode(RNode *neighb, RNode **edpp, RNode *newn, int level, float *vert)
404 {
405         /*      IF neighb->ed points to newn->up
406          *              IF edgelevels equal
407                                 IF testvert is in neighb->ed
408                                         change pointers both ways
409                                 ELSE
410                                         RETURN
411                         ELSE
412                                 IF neighb edgelevel is deeper
413                                         change neighb pointer
414                 
415          */
416         int ok= 0;
417         
418         if(neighb==0) return;
419         if(newn->up==0) return;
420         
421         if(neighb->ed1==newn->up) {
422                 if(neighb->lev1==level) {
423                         if(vert==neighb->v1 || vert==neighb->v2) {
424                                 *edpp= neighb;
425                                 neighb->ed1= newn;
426                         }
427                         else return;
428                 }
429                 else if(neighb->lev1>level) {
430                         neighb->ed1= newn;
431                 }
432                 ok= 1;
433         }
434         else if(neighb->ed2==newn->up) {
435                 if(neighb->lev2==level) {
436                         if(vert==neighb->v2 || vert==neighb->v3) {
437                                 *edpp= neighb;
438                                 neighb->ed2= newn;
439                         }
440                         else return;
441                 }
442                 else if(neighb->lev2>level) {
443                         neighb->ed2= newn;
444                 }
445                 ok= 1;
446         }
447         else if(neighb->ed3==newn->up) {
448                 if(neighb->lev3==level) {
449                         if(neighb->type==3) {
450                                 if(vert==neighb->v3 || vert==neighb->v1) {
451                                         *edpp= neighb;
452                                         neighb->ed3= newn;
453                                 }
454                                 else return;
455                         }
456                         else {
457                                 if(vert==neighb->v3 || vert==neighb->v4) {
458                                         *edpp= neighb;
459                                         neighb->ed3= newn;
460                                 }
461                                 else return;
462                         }
463                 }
464                 else if(neighb->lev3>level) {
465                         neighb->ed3= newn;
466                 }
467                 ok= 1;
468         }
469         else if(neighb->ed4==newn->up) {
470                 if(neighb->lev4==level) {
471                         if(vert==neighb->v4 || vert==neighb->v1) {
472                                 *edpp= neighb;
473                                 neighb->ed4= newn;
474                         }
475                         else return;
476                 }
477                 else if(neighb->lev4>level) {
478                         neighb->ed4= newn;
479                 }
480                 ok= 1;
481         }
482         
483         if(ok && neighb->down1) {
484                 replaceTestNode(neighb->down1, edpp, newn, level, vert);
485                 replaceTestNode(neighb->down2, edpp, newn, level, vert);
486         }
487         
488 }
489
490 int setvertexpointersNode(RNode *neighb, RNode *node, int level, float **v1, float **v2)
491 {
492         /* compares edgelevels , if equal it sets the vertexpointers */
493         
494         if(neighb==0) return 0;
495         
496         if(neighb->ed1==node) {
497                 if(neighb->lev1==level) {
498                         *v1= neighb->v1;
499                         *v2= neighb->v2;
500                         return 1;
501                 }
502         }
503         else if(neighb->ed2==node) {
504                 if(neighb->lev2==level) {
505                         *v1= neighb->v2;
506                         *v2= neighb->v3;
507                         return 1;
508                 }
509         }
510         else if(neighb->ed3==node) {
511                 if(neighb->lev3==level) {
512                         if(neighb->type==3) {
513                                 *v1= neighb->v3;
514                                 *v2= neighb->v1;
515                         }
516                         else {
517                                 *v1= neighb->v3;
518                                 *v2= neighb->v4;
519                         }
520                         return 1;
521                 }
522         }
523         else if(neighb->ed4==node) {
524                 if(neighb->lev4==level) {
525                         *v1= neighb->v4;
526                         *v2= neighb->v1;
527                         return 1;
528                 }
529         }
530         return 0;
531 }
532
533 float edlen(float *v1, float *v2)
534 {
535         return (v1[0]-v2[0])*(v1[0]-v2[0])+ (v1[1]-v2[1])*(v1[1]-v2[1])+ (v1[2]-v2[2])*(v1[2]-v2[2]);
536 }
537
538
539 void subdivideTriNode(RNode *node, RNode *edge)
540 {
541         RNode *n1, *n2, *up;
542         float fu, fv, fl, *v1, *v2; /*  , AreaT3Dfl(); ... from arithb... */
543         int uvl;
544         
545         if(node->down1 || node->down2) {
546                 /* printf("trinode: subd already done\n"); */
547                 return;
548         }
549         
550         /* defines subdivide direction */
551
552         if(edge==0) {
553                 /* areathreshold */
554                 if(node->area<nodelimit) return;
555
556                 fu= edlen(node->v1, node->v2);
557                 fv= edlen(node->v2, node->v3);
558                 fl= edlen(node->v3, node->v1);
559
560                 if(fu>fv && fu>fl) uvl= 1;
561                 else if(fv>fu && fv>fl) uvl= 2;
562                 else uvl= 3;
563         }
564         else {
565                 
566                 if(edge==node->ed1) uvl= 1;
567                 else if(edge==node->ed2) uvl= 2;
568                 else uvl= 3;
569         }
570         
571         /*  should neighbour nodes be deeper? Recursive! */
572         n1= 0; 
573         if(uvl==1) {
574                 if(node->ed1 && node->ed1->down1==0) n1= node->ed1;
575         }
576         else if(uvl==2) {
577                 if(node->ed2 && node->ed2->down1==0) n1= node->ed2;
578         }
579         else {
580                 if(node->ed3 && node->ed3->down1==0) n1= node->ed3;
581         }
582         if(n1) {
583                 up= node->up;
584                 while(up) {                                                     /* also test for ed4 !!! */
585                         if(n1->ed1==up || n1->ed2==up || n1->ed3==up || n1->ed4==up) {
586                                 subdivideNode(n1, up);
587                                 break;
588                         }
589                         up= up->up;
590                 }
591         }
592         
593         /* the subdividing */
594         n1= mallocNode();
595         memcpy(n1, node, sizeof(RNode));
596         n2= mallocNode();
597         memcpy(n2, node, sizeof(RNode));
598
599         n1->up= node;
600         n2->up= node;
601         
602         node->down1= n1;
603         node->down2= n2;
604
605         /* subdivide edge 1 */
606         if(uvl==1) {
607         
608                 /* FIRST NODE  gets edge 2 */
609                 n1->ed3= n2;
610                 n1->lev3= 0;
611                 replaceAllNode(n1->ed2, n1);
612                 n1->lev1++;
613                 replaceTestNode(n1->ed1, &(n1->ed1), n1, n1->lev1, n1->v2);
614
615                 /* SECOND NODE  gets edge 3 */
616                 n2->ed2= n1;
617                 n2->lev2= 0;
618                 replaceAllNode(n2->ed3, n2);
619                 n2->lev1++;
620                 replaceTestNode(n2->ed1, &(n2->ed1), n2, n2->lev1, n2->v1);
621                 
622                 /* NEW VERTEX from edge 1 */
623                 if( setvertexpointersNode(n1->ed1, n1, n1->lev1, &v1, &v2) ) {  /* nodes have equal levels */
624                         if(v1== n1->v2) {
625                                 n1->v1= v2;
626                                 n2->v2= v2;
627                         }
628                         else {
629                                 n1->v1= v1;
630                                 n2->v2= v1;
631                         }
632                 }
633                 else {
634                         n1->v1= n2->v2= mallocVert();
635                         n1->v1[0]= 0.5*(node->v1[0]+ node->v2[0]);
636                         n1->v1[1]= 0.5*(node->v1[1]+ node->v2[1]);
637                         n1->v1[2]= 0.5*(node->v1[2]+ node->v2[2]);
638                         n1->v1[3]= node->v1[3]; /* color */
639                 }
640         }
641         else if(uvl==2) {
642         
643                 /* FIRST NODE gets edge 1 */
644                 n1->ed3= n2;
645                 n1->lev3= 0;
646                 replaceAllNode(n1->ed1, n1);
647                 n1->lev2++;
648                 replaceTestNode(n1->ed2, &(n1->ed2), n1, n1->lev2, n1->v2);
649
650                 /* SECOND NODE gets edge 3 */
651                 n2->ed1= n1;
652                 n2->lev1= 0;
653                 replaceAllNode(n2->ed3, n2);
654                 n2->lev2++;
655                 replaceTestNode(n2->ed2, &(n2->ed2), n2, n2->lev2, n2->v3);
656                 
657                 /* NEW VERTEX from edge 2 */
658                 if( setvertexpointersNode(n1->ed2, n1, n1->lev2, &v1, &v2) ) {  /* nodes have equal levels */
659                         if(v1== n1->v2) {
660                                 n1->v3= v2;
661                                 n2->v2= v2;
662                         }
663                         else {
664                                 n1->v3= v1;
665                                 n2->v2= v1;
666                         }
667                 }
668                 else {
669                         n1->v3= n2->v2= mallocVert();
670                         n1->v3[0]= 0.5*(node->v2[0]+ node->v3[0]);
671                         n1->v3[1]= 0.5*(node->v2[1]+ node->v3[1]);
672                         n1->v3[2]= 0.5*(node->v2[2]+ node->v3[2]);
673                         n1->v3[3]= node->v1[3]; /* color */
674                 }
675         }
676         else if(uvl==3) {
677         
678                 /* FIRST NODE gets edge 1 */
679                 n1->ed2= n2;
680                 n1->lev2= 0;
681                 replaceAllNode(n1->ed1, n1);
682                 n1->lev3++;
683                 replaceTestNode(n1->ed3, &(n1->ed3), n1, n1->lev3, n1->v1);
684
685                 /* SECOND NODE gets edge 2 */
686                 n2->ed1= n1;
687                 n2->lev1= 0;
688                 replaceAllNode(n2->ed2, n2);
689                 n2->lev3++;
690                 replaceTestNode(n2->ed3, &(n2->ed3), n2, n2->lev3, n2->v3);
691                 
692                 /* NEW VERTEX from edge 3 */
693                 if( setvertexpointersNode(n1->ed3, n1, n1->lev3, &v1, &v2) ) {  /* nodes have equal levels */
694                         if(v1== n1->v1) {
695                                 n1->v3= v2;
696                                 n2->v1= v2;
697                         }
698                         else {
699                                 n1->v3= v1;
700                                 n2->v1= v1;
701                         }
702                 }
703                 else {
704                         n1->v3= n2->v1= mallocVert();
705                         n1->v3[0]= 0.5*(node->v1[0]+ node->v3[0]);
706                         n1->v3[1]= 0.5*(node->v1[1]+ node->v3[1]);
707                         n1->v3[2]= 0.5*(node->v1[2]+ node->v3[2]);
708                         n1->v3[3]= node->v3[3]; /* color */
709                 }
710         }
711         n1->area= AreaT3Dfl(n1->v1, n1->v2, n1->v3);
712         n2->area= AreaT3Dfl(n2->v1, n2->v2, n2->v3);
713
714 }
715
716
717 void subdivideNode(RNode *node, RNode *edge)
718 {
719         RNode *n1, *n2, *up;
720         float fu, fv, *v1, *v2;/*, AreaQ3Dfl(); ... from arithb... */
721         int uvl;
722         
723         if(Ntotnode>RG.maxnode) return;
724         
725         if(node->type==3) {
726                 subdivideTriNode(node, edge);
727                 return;
728         }
729
730         if(node->down1 || node->down2) {
731                 /* printf("subdivide Node: already done \n"); */
732                 return;
733         }
734         
735         /* defines subdivide direction */
736
737         if(edge==0) {
738                 /* areathreshold */
739                 if(node->area<nodelimit) {
740                         return;
741                 }
742                 fu= fabs(node->v1[0]- node->v2[0])+ fabs(node->v1[1]- node->v2[1]) +fabs(node->v1[2]- node->v2[2]);
743                 fv= fabs(node->v1[0]- node->v4[0])+ fabs(node->v1[1]- node->v4[1]) +fabs(node->v1[2]- node->v4[2]);
744                 if(fu>fv) uvl= 1;
745                 else uvl= 2;
746         }
747         else {
748                 if(edge==node->ed1 || edge==node->ed3) uvl= 1;
749                 else uvl= 2;
750         }
751         
752         /*  do neighbour nodes have to be deeper? Recursive! */
753         n1= n2= 0; 
754         if(uvl==1) {
755                 if(node->ed1 && node->ed1->down1==0) n1= node->ed1;
756                 if(node->ed3 && node->ed3->down1==0) n2= node->ed3;
757         }
758         else {
759                 if(node->ed2 && node->ed2->down1==0) n1= node->ed2;
760                 if(node->ed4 && node->ed4->down1==0) n2= node->ed4;
761         }
762         if(n1) {
763                 up= node->up;
764                 while(up) {
765                         if(n1->ed1==up || n1->ed2==up || n1->ed3==up || n1->ed4==up) {
766                                 /* printf("recurs subd\n"); */
767                                 subdivideNode(n1, up);
768                                 break;
769                         }
770                         up= up->up;
771                 }
772         }
773         if(n2) {
774                 up= node->up;
775                 while(up) {
776                         if(n2->ed1==up || n2->ed2==up || n2->ed3==up || n2->ed4==up) {
777                                 /* printf("recurs subd\n"); */
778                                 subdivideNode(n2, up);
779                                 break;
780                         }
781                         up= up->up;
782                 }
783         }
784
785         /* the subdividing */
786         n1= mallocNode();
787         memcpy(n1, node, sizeof(RNode));
788         n2= mallocNode();
789         memcpy(n2, node, sizeof(RNode));
790
791         n1->up= node;
792         n2->up= node;
793         
794         node->down1= n1;
795         node->down2= n2;
796
797         /* subdivide edge 1 and 3 */
798         if(uvl==1) {
799                 
800                 /* FIRST NODE  gets edge 2 */
801                 n1->ed4= n2;
802                 n1->lev4= 0;
803                 replaceAllNode(n1->ed2, n1);
804                 n1->lev1++;
805                 n1->lev3++;
806                 replaceTestNode(n1->ed1, &(n1->ed1), n1, n1->lev1, n1->v2);
807                 replaceTestNode(n1->ed3, &(n1->ed3), n1, n1->lev3, n1->v3);
808
809                 /* SECOND NODE gets edge 4 */
810                 n2->ed2= n1;
811                 n2->lev2= 0;
812                 replaceAllNode(n2->ed4, n2);
813                 n2->lev1++;
814                 n2->lev3++;
815                 replaceTestNode(n2->ed1, &(n2->ed1), n2, n2->lev1, n2->v1);
816                 replaceTestNode(n2->ed3, &(n2->ed3), n2, n2->lev3, n2->v4);
817                 
818                 /* NEW VERTEX from edge 1 */
819                 if( setvertexpointersNode(n1->ed1, n1, n1->lev1, &v1, &v2) ) {  /* nodes have equal levels */
820                         if(v1== n1->v2) {
821                                 n1->v1= v2;
822                                 n2->v2= v2;
823                         }
824                         else {
825                                 n1->v1= v1;
826                                 n2->v2= v1;
827                         }
828                 }
829                 else {
830                         n1->v1= n2->v2= mallocVert();
831                         n1->v1[0]= 0.5*(node->v1[0]+ node->v2[0]);
832                         n1->v1[1]= 0.5*(node->v1[1]+ node->v2[1]);
833                         n1->v1[2]= 0.5*(node->v1[2]+ node->v2[2]);
834                         n1->v1[3]= node->v1[3]; /* color */
835                 }
836                 
837                 /* NEW VERTEX from edge 3 */
838                 if( setvertexpointersNode(n1->ed3, n1, n1->lev3, &v1, &v2) ) {  /* nodes have equal levels */
839                         if(v1== n1->v3) {
840                                 n1->v4= v2;
841                                 n2->v3= v2;
842                         }
843                         else {
844                                 n1->v4= v1;
845                                 n2->v3= v1;
846                         }
847                 }
848                 else {
849                         n1->v4= n2->v3= mallocVert();
850                         n1->v4[0]= 0.5*(node->v3[0]+ node->v4[0]);
851                         n1->v4[1]= 0.5*(node->v3[1]+ node->v4[1]);
852                         n1->v4[2]= 0.5*(node->v3[2]+ node->v4[2]);
853                         n1->v4[3]= node->v4[3]; /* color */
854                 }
855         }
856         /* subdivide edge 2 and 4 */
857         else if(uvl==2) {
858                 
859                 /* FIRST NODE gets edge 1 */
860                 n1->ed3= n2;
861                 n1->lev3= 0;
862                 replaceAllNode(n1->ed1, n1);
863                 n1->lev2++;
864                 n1->lev4++;
865                 replaceTestNode(n1->ed2, &(n1->ed2), n1, n1->lev2, n1->v2);
866                 replaceTestNode(n1->ed4, &(n1->ed4), n1, n1->lev4, n1->v1);
867
868                 /* SECOND NODE gets edge 3 */
869                 n2->ed1= n1;
870                 n2->lev1= 0;
871                 replaceAllNode(n2->ed3, n2);
872                 n2->lev2++;
873                 n2->lev4++;
874                 replaceTestNode(n2->ed2, &(n2->ed2), n2, n2->lev2, n2->v3);
875                 replaceTestNode(n2->ed4, &(n2->ed4), n2, n2->lev4, n2->v4);
876
877                 /* NEW VERTEX from edge 2 */
878                 if( setvertexpointersNode(n1->ed2, n1, n1->lev2, &v1, &v2) ) {  /* nodes have equal levels */
879                         if(v1== n1->v2) {
880                                 n1->v3= v2;
881                                 n2->v2= v2;
882                         }
883                         else {
884                                 n1->v3= v1;
885                                 n2->v2= v1;
886                         }
887                 }
888                 else {
889                         n1->v3= n2->v2= mallocVert();
890                         n1->v3[0]= 0.5*(node->v2[0]+ node->v3[0]);
891                         n1->v3[1]= 0.5*(node->v2[1]+ node->v3[1]);
892                         n1->v3[2]= 0.5*(node->v2[2]+ node->v3[2]);
893                         n1->v3[3]= node->v3[3]; /* color */
894                 }
895
896                 /* NEW VERTEX from edge 4 */
897                 if( setvertexpointersNode(n1->ed4, n1, n1->lev4, &v1, &v2) ) {  /* nodes have equal levels */
898                         if(v1== n1->v1) {
899                                 n1->v4= v2;
900                                 n2->v1= v2;
901                         }
902                         else {
903                                 n1->v4= v1;
904                                 n2->v1= v1;
905                         }
906                 }
907                 else {
908                         n1->v4= n2->v1= mallocVert();
909                         n1->v4[0]= 0.5*(node->v1[0]+ node->v4[0]);
910                         n1->v4[1]= 0.5*(node->v1[1]+ node->v4[1]);
911                         n1->v4[2]= 0.5*(node->v1[2]+ node->v4[2]);
912                         n1->v4[3]= node->v4[3]; /* color */
913                 }
914         }
915         
916         n1->area= AreaQ3Dfl(n1->v1, n1->v2, n1->v3, n1->v4);
917         n2->area= AreaQ3Dfl(n2->v1, n2->v2, n2->v3, n2->v4);
918
919 }
920
921 int comparelevel(RNode *node, RNode *nb, int level)
922 {
923         /* recursive descent: test with deepest node */
924         /* return 1 means equal or higher */
925         
926         if(nb==0) return 1;
927         
928         if(nb->down1) {
929                 return 0;
930                 
931                 /*              THERE IS AN ERROR HERE, BUT WHAT?  (without this function the system
932                         works too, but is slower) (ton) */
933
934                 /*
935                 n1= nb->down1;
936                 if(n1->ed1==node) return comparelevel(node, n1, level);
937                 if(n1->ed2==node) return comparelevel(node, n1, level);
938                 if(n1->ed3==node) return comparelevel(node, n1, level);
939                 if(n1->ed4==node) return comparelevel(node, n1, level);
940                 n1= nb->down2;
941                 if(n1->ed1==node) return comparelevel(node, n1, level);
942                 if(n1->ed2==node) return comparelevel(node, n1, level);
943                 if(n1->ed3==node) return comparelevel(node, n1, level);
944                 if(n1->ed4==node) return comparelevel(node, n1, level);
945                 printf(" dit kan niet ");
946                 return 0;
947                 */
948                 
949         }
950         
951         if(nb->down1==0) {
952                 /* if(nb->ed1==node) return (nb->lev1<=level); */
953                 /* if(nb->ed2==node) return (nb->lev2<=level); */
954                 /* if(nb->ed3==node) return (nb->lev3<=level); */
955                 /* if(nb->ed4==node) return (nb->lev4<=level); */
956                 
957                 return 1;       /* is higher node */
958         }
959         return 1;
960 }
961
962 void deleteTriNodes(RNode *node)        /* both children of node */
963 {
964         RNode *n1, *n2;
965         
966         /* if neighbour nodes are deeper: no delete */
967         /* just test 2 nodes, from the others the level doesn't change */
968         
969         n1= node->down1;
970         n2= node->down2;
971
972         if(n1==0 || n2==0) return;
973         
974         if(n1->down1 || n2->down1) return;
975         
976         /* at the edges no subdivided node is allowed */
977
978         if(n1->ed1 && n1->ed1->down1) return;
979         if(n1->ed2 && n1->ed2->down1) return;
980         if(n1->ed3 && n1->ed3->down1) return;
981
982         if(n2->ed1 && n2->ed1->down1) return;
983         if(n2->ed2 && n2->ed2->down1) return;
984         if(n2->ed3 && n2->ed3->down1) return;
985                                 
986         replaceAllNodeInv(n1->ed1, n1);
987         replaceAllNodeInv(n1->ed2, n1);
988         replaceAllNodeInv(n1->ed3, n1);
989
990         replaceAllNodeUp(n1->ed1, n1);
991         replaceAllNodeUp(n1->ed2, n1);
992         replaceAllNodeUp(n1->ed3, n1);
993
994         replaceAllNodeInv(n2->ed1, n2);
995         replaceAllNodeInv(n2->ed2, n2);
996         replaceAllNodeInv(n2->ed3, n2);
997
998         replaceAllNodeUp(n2->ed1, n2);
999         replaceAllNodeUp(n2->ed2, n2);
1000         replaceAllNodeUp(n2->ed3, n2);
1001
1002         n1->down1= (RNode *)12; /* for debug */
1003         n2->down1= (RNode *)12;
1004         
1005         freeNode(n1);
1006         freeNode(n2);
1007         node->down1= node->down2= 0;
1008         
1009 }
1010
1011         /* both children of node */
1012 void deleteNodes(RNode *node)
1013 {
1014         RNode *n1, *n2;
1015
1016         /* if neighbour nodes are deeper: no delete */
1017         /* just test 2 nodes, from the others the level doesn't change */
1018
1019         if(node->type==3) {
1020                 deleteTriNodes(node);
1021                 return;
1022         }
1023         
1024         n1= node->down1;
1025         n2= node->down2;
1026
1027         if(n1==0 || n2==0) return;
1028         
1029         if(n1->down1 || n2->down1) return;
1030         
1031         if(n1->ed3==n2) {
1032
1033                 /* at the edges no subdivided node is allowed */
1034
1035                 if(n1->ed1 && n1->ed1->down1) return;
1036                 if(n1->ed2 && n1->ed2->down1) return;
1037                 if(n1->ed4 && n1->ed4->down1) return;
1038
1039                 if(n2->ed2 && n2->ed2->down1) return;
1040                 if(n2->ed3 && n2->ed3->down1) return;
1041                 if(n2->ed4 && n2->ed4->down1) return;
1042                                         
1043                 replaceAllNodeInv(n1->ed1, n1);
1044                 replaceAllNodeInv(n1->ed2, n1);
1045                 replaceAllNodeInv(n1->ed4, n1);
1046
1047                 replaceAllNodeUp(n1->ed1, n1);
1048                 replaceAllNodeUp(n1->ed2, n1);
1049                 replaceAllNodeUp(n1->ed4, n1);
1050
1051                 replaceAllNodeInv(n2->ed2, n2);
1052                 replaceAllNodeInv(n2->ed3, n2);
1053                 replaceAllNodeInv(n2->ed4, n2);
1054
1055                 replaceAllNodeUp(n2->ed2, n2);
1056                 replaceAllNodeUp(n2->ed3, n2);
1057                 replaceAllNodeUp(n2->ed4, n2);
1058
1059                 n1->down1= (RNode *)12; /* for debug */
1060                 n2->down1= (RNode *)12;
1061                 
1062                 freeNode(n1);
1063                 freeNode(n2);
1064                 node->down1= node->down2= 0;
1065                 
1066                 return;
1067         }
1068         else if(n1->ed4==n2) {
1069
1070                 if(n1->ed1 && n1->ed1->down1) return;
1071                 if(n1->ed2 && n1->ed2->down1) return;
1072                 if(n1->ed3 && n1->ed3->down1) return;
1073
1074                 if(n2->ed1 && n2->ed1->down1) return;
1075                 if(n2->ed3 && n2->ed3->down1) return;
1076                 if(n2->ed4 && n2->ed4->down1) return;
1077                                         
1078                 replaceAllNodeInv(n1->ed1, n1);
1079                 replaceAllNodeInv(n1->ed2, n1);
1080                 replaceAllNodeInv(n1->ed3, n1);
1081
1082                 replaceAllNodeUp(n1->ed1, n1);
1083                 replaceAllNodeUp(n1->ed2, n1);
1084                 replaceAllNodeUp(n1->ed3, n1);
1085
1086                 replaceAllNodeInv(n2->ed1, n2);
1087                 replaceAllNodeInv(n2->ed3, n2);
1088                 replaceAllNodeInv(n2->ed4, n2);
1089
1090                 replaceAllNodeUp(n2->ed1, n2);
1091                 replaceAllNodeUp(n2->ed3, n2);
1092                 replaceAllNodeUp(n2->ed4, n2);
1093
1094                 n1->down1= (RNode *)12; /* for debug */
1095                 n2->down1= (RNode *)12;
1096                 
1097                 freeNode(n1);
1098                 freeNode(n2);
1099                 node->down1= node->down2= 0;
1100                 
1101                 return;
1102         }
1103
1104 }
1105
1106