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