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