15901a66f1ed7f5e4d5557eb793b246cc4894f3f
[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 RNode *mallocNode()
263 {
264         Ntotnode++;
265         return (RNode *)malloc_fast(sizeof(RNode));
266 }
267
268 RNode *callocNode()
269 {
270         Ntotnode++;
271         return (RNode *)calloc_fast(sizeof(RNode));
272 }
273
274 void freeNode(RNode *node)
275 {
276         free_fast(node, sizeof(RNode));
277         Ntotnode--;
278 }
279
280 void freeNode_recurs(RNode *node)
281 {
282
283         if(node->down1) {
284                 freeNode_recurs(node->down1);
285                 freeNode_recurs(node->down2);
286         }
287
288         node->down1= node->down2= 0;
289         freeNode(node);
290
291 }
292
293 RPatch *mallocPatch()
294 {
295         Ntotpatch++;
296         return (RPatch *)malloc_fast(sizeof(RPatch));
297 }
298
299 RPatch *callocPatch()
300 {
301         Ntotpatch++;
302         return (RPatch *)calloc_fast(sizeof(RPatch));
303 }
304
305 void freePatch(RPatch *patch)
306 {
307         free_fast(patch, sizeof(RPatch));
308         Ntotpatch--;
309 }
310
311 /* ************  SUBDIVIDE  ***********  */
312
313
314 void replaceAllNode(RNode *neighb, RNode *newn)
315 {
316         /* changes from all neighbours the edgepointers that point to newn->up in new */
317         int ok= 0;
318         
319         
320         if(neighb==0) return;
321         if(newn->up==0) return;
322         
323         if(neighb->ed1==newn->up) {
324                 neighb->ed1= newn;
325                 ok= 1;
326         }
327         else if(neighb->ed2==newn->up) {
328                 neighb->ed2= newn;
329                 ok= 1;
330         }
331         else if(neighb->ed3==newn->up) {
332                 neighb->ed3= newn;
333                 ok= 1;
334         }
335         else if(neighb->ed4==newn->up) {
336                 neighb->ed4= newn;
337                 ok= 1;
338         }
339         
340         if(ok && neighb->down1) {
341                 replaceAllNode(neighb->down1, newn);
342                 replaceAllNode(neighb->down2, newn);
343         }
344 }
345
346 void replaceAllNodeInv(RNode *neighb, RNode *old)
347 {
348         /* changes from all neighbours the edgepointers that point to old in old->up */
349         if(neighb==0) return;
350         if(old->up==0) return;
351         
352         if(neighb->ed1==old) {
353                 neighb->ed1= old->up;
354         }
355         else if(neighb->ed2==old) {
356                 neighb->ed2= old->up;
357         }
358         else if(neighb->ed3==old) {
359                 neighb->ed3= old->up;
360         }
361         else if(neighb->ed4==old) {
362                 neighb->ed4= old->up;
363         }
364         
365         if(neighb->down1) {
366                 replaceAllNodeInv(neighb->down1, old);
367                 replaceAllNodeInv(neighb->down2, old);
368         }
369 }
370
371 void replaceAllNodeUp(RNode *neighb, RNode *old)
372 {
373         /* changes from all neighbours the edgepointers that point to old in old->up */
374         if(neighb==0) return;
375         if(old->up==0) return;
376         neighb= neighb->up;
377         if(neighb==0) return;
378         
379         if(neighb->ed1==old) {
380                 neighb->ed1= old->up;
381         }
382         else if(neighb->ed2==old) {
383                 neighb->ed2= old->up;
384         }
385         else if(neighb->ed3==old) {
386                 neighb->ed3= old->up;
387         }
388         else if(neighb->ed4==old) {
389                 neighb->ed4= old->up;
390         }
391         
392         if(neighb->up) {
393                 replaceAllNodeUp(neighb, old);
394         }
395 }
396
397
398 void replaceTestNode(RNode *neighb, RNode **edpp, RNode *newn, int level, float *vert)
399 {
400         /*      IF neighb->ed points to newn->up
401          *              IF edgelevels equal
402                                 IF testvert is in neighb->ed
403                                         change pointers both ways
404                                 ELSE
405                                         RETURN
406                         ELSE
407                                 IF neighb edgelevel is deeper
408                                         change neighb pointer
409                 
410          */
411         int ok= 0;
412         
413         if(neighb==0) return;
414         if(newn->up==0) return;
415         
416         if(neighb->ed1==newn->up) {
417                 if(neighb->lev1==level) {
418                         if(vert==neighb->v1 || vert==neighb->v2) {
419                                 *edpp= neighb;
420                                 neighb->ed1= newn;
421                         }
422                         else return;
423                 }
424                 else if(neighb->lev1>level) {
425                         neighb->ed1= newn;
426                 }
427                 ok= 1;
428         }
429         else if(neighb->ed2==newn->up) {
430                 if(neighb->lev2==level) {
431                         if(vert==neighb->v2 || vert==neighb->v3) {
432                                 *edpp= neighb;
433                                 neighb->ed2= newn;
434                         }
435                         else return;
436                 }
437                 else if(neighb->lev2>level) {
438                         neighb->ed2= newn;
439                 }
440                 ok= 1;
441         }
442         else if(neighb->ed3==newn->up) {
443                 if(neighb->lev3==level) {
444                         if(neighb->type==3) {
445                                 if(vert==neighb->v3 || vert==neighb->v1) {
446                                         *edpp= neighb;
447                                         neighb->ed3= newn;
448                                 }
449                                 else return;
450                         }
451                         else {
452                                 if(vert==neighb->v3 || vert==neighb->v4) {
453                                         *edpp= neighb;
454                                         neighb->ed3= newn;
455                                 }
456                                 else return;
457                         }
458                 }
459                 else if(neighb->lev3>level) {
460                         neighb->ed3= newn;
461                 }
462                 ok= 1;
463         }
464         else if(neighb->ed4==newn->up) {
465                 if(neighb->lev4==level) {
466                         if(vert==neighb->v4 || vert==neighb->v1) {
467                                 *edpp= neighb;
468                                 neighb->ed4= newn;
469                         }
470                         else return;
471                 }
472                 else if(neighb->lev4>level) {
473                         neighb->ed4= newn;
474                 }
475                 ok= 1;
476         }
477         
478         if(ok && neighb->down1) {
479                 replaceTestNode(neighb->down1, edpp, newn, level, vert);
480                 replaceTestNode(neighb->down2, edpp, newn, level, vert);
481         }
482         
483 }
484
485 int setvertexpointersNode(RNode *neighb, RNode *node, int level, float **v1, float **v2)
486 {
487         /* compares edgelevels , if equal it sets the vertexpointers */
488         
489         if(neighb==0) return 0;
490         
491         if(neighb->ed1==node) {
492                 if(neighb->lev1==level) {
493                         *v1= neighb->v1;
494                         *v2= neighb->v2;
495                         return 1;
496                 }
497         }
498         else if(neighb->ed2==node) {
499                 if(neighb->lev2==level) {
500                         *v1= neighb->v2;
501                         *v2= neighb->v3;
502                         return 1;
503                 }
504         }
505         else if(neighb->ed3==node) {
506                 if(neighb->lev3==level) {
507                         if(neighb->type==3) {
508                                 *v1= neighb->v3;
509                                 *v2= neighb->v1;
510                         }
511                         else {
512                                 *v1= neighb->v3;
513                                 *v2= neighb->v4;
514                         }
515                         return 1;
516                 }
517         }
518         else if(neighb->ed4==node) {
519                 if(neighb->lev4==level) {
520                         *v1= neighb->v4;
521                         *v2= neighb->v1;
522                         return 1;
523                 }
524         }
525         return 0;
526 }
527
528 float edlen(float *v1, float *v2)
529 {
530         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]);
531 }
532
533
534 void subdivideTriNode(RNode *node, RNode *edge)
535 {
536         RNode *n1, *n2, *up;
537         float fu, fv, fl, *v1, *v2; /*  , AreaT3Dfl(); ... from arithb... */
538         int uvl;
539         
540         if(node->down1 || node->down2) {
541                 /* printf("trinode: subd already done\n"); */
542                 return;
543         }
544         
545         /* defines subdivide direction */
546
547         if(edge==0) {
548                 /* areathreshold */
549                 if(node->area<nodelimit) return;
550
551                 fu= edlen(node->v1, node->v2);
552                 fv= edlen(node->v2, node->v3);
553                 fl= edlen(node->v3, node->v1);
554
555                 if(fu>fv && fu>fl) uvl= 1;
556                 else if(fv>fu && fv>fl) uvl= 2;
557                 else uvl= 3;
558         }
559         else {
560                 
561                 if(edge==node->ed1) uvl= 1;
562                 else if(edge==node->ed2) uvl= 2;
563                 else uvl= 3;
564         }
565         
566         /*  should neighbour nodes be deeper? Recursive! */
567         n1= 0; 
568         if(uvl==1) {
569                 if(node->ed1 && node->ed1->down1==0) n1= node->ed1;
570         }
571         else if(uvl==2) {
572                 if(node->ed2 && node->ed2->down1==0) n1= node->ed2;
573         }
574         else {
575                 if(node->ed3 && node->ed3->down1==0) n1= node->ed3;
576         }
577         if(n1) {
578                 up= node->up;
579                 while(up) {                                                     /* also test for ed4 !!! */
580                         if(n1->ed1==up || n1->ed2==up || n1->ed3==up || n1->ed4==up) {
581                                 subdivideNode(n1, up);
582                                 break;
583                         }
584                         up= up->up;
585                 }
586         }
587         
588         /* the subdividing */
589         n1= mallocNode();
590         memcpy(n1, node, sizeof(RNode));
591         n2= mallocNode();
592         memcpy(n2, node, sizeof(RNode));
593
594         n1->up= node;
595         n2->up= node;
596         
597         node->down1= n1;
598         node->down2= n2;
599
600         /* subdivide edge 1 */
601         if(uvl==1) {
602         
603                 /* FIRST NODE  gets edge 2 */
604                 n1->ed3= n2;
605                 n1->lev3= 0;
606                 replaceAllNode(n1->ed2, n1);
607                 n1->lev1++;
608                 replaceTestNode(n1->ed1, &(n1->ed1), n1, n1->lev1, n1->v2);
609
610                 /* SECOND NODE  gets edge 3 */
611                 n2->ed2= n1;
612                 n2->lev2= 0;
613                 replaceAllNode(n2->ed3, n2);
614                 n2->lev1++;
615                 replaceTestNode(n2->ed1, &(n2->ed1), n2, n2->lev1, n2->v1);
616                 
617                 /* NEW VERTEX from edge 1 */
618                 if( setvertexpointersNode(n1->ed1, n1, n1->lev1, &v1, &v2) ) {  /* nodes have equal levels */
619                         if(v1== n1->v2) {
620                                 n1->v1= v2;
621                                 n2->v2= v2;
622                         }
623                         else {
624                                 n1->v1= v1;
625                                 n2->v2= v1;
626                         }
627                 }
628                 else {
629                         n1->v1= n2->v2= mallocVert();
630                         n1->v1[0]= 0.5*(node->v1[0]+ node->v2[0]);
631                         n1->v1[1]= 0.5*(node->v1[1]+ node->v2[1]);
632                         n1->v1[2]= 0.5*(node->v1[2]+ node->v2[2]);
633                         n1->v1[3]= node->v1[3]; /* color */
634                 }
635         }
636         else if(uvl==2) {
637         
638                 /* FIRST NODE gets edge 1 */
639                 n1->ed3= n2;
640                 n1->lev3= 0;
641                 replaceAllNode(n1->ed1, n1);
642                 n1->lev2++;
643                 replaceTestNode(n1->ed2, &(n1->ed2), n1, n1->lev2, n1->v2);
644
645                 /* SECOND NODE gets edge 3 */
646                 n2->ed1= n1;
647                 n2->lev1= 0;
648                 replaceAllNode(n2->ed3, n2);
649                 n2->lev2++;
650                 replaceTestNode(n2->ed2, &(n2->ed2), n2, n2->lev2, n2->v3);
651                 
652                 /* NEW VERTEX from edge 2 */
653                 if( setvertexpointersNode(n1->ed2, n1, n1->lev2, &v1, &v2) ) {  /* nodes have equal levels */
654                         if(v1== n1->v2) {
655                                 n1->v3= v2;
656                                 n2->v2= v2;
657                         }
658                         else {
659                                 n1->v3= v1;
660                                 n2->v2= v1;
661                         }
662                 }
663                 else {
664                         n1->v3= n2->v2= mallocVert();
665                         n1->v3[0]= 0.5*(node->v2[0]+ node->v3[0]);
666                         n1->v3[1]= 0.5*(node->v2[1]+ node->v3[1]);
667                         n1->v3[2]= 0.5*(node->v2[2]+ node->v3[2]);
668                         n1->v3[3]= node->v1[3]; /* color */
669                 }
670         }
671         else if(uvl==3) {
672         
673                 /* FIRST NODE gets edge 1 */
674                 n1->ed2= n2;
675                 n1->lev2= 0;
676                 replaceAllNode(n1->ed1, n1);
677                 n1->lev3++;
678                 replaceTestNode(n1->ed3, &(n1->ed3), n1, n1->lev3, n1->v1);
679
680                 /* SECOND NODE gets edge 2 */
681                 n2->ed1= n1;
682                 n2->lev1= 0;
683                 replaceAllNode(n2->ed2, n2);
684                 n2->lev3++;
685                 replaceTestNode(n2->ed3, &(n2->ed3), n2, n2->lev3, n2->v3);
686                 
687                 /* NEW VERTEX from edge 3 */
688                 if( setvertexpointersNode(n1->ed3, n1, n1->lev3, &v1, &v2) ) {  /* nodes have equal levels */
689                         if(v1== n1->v1) {
690                                 n1->v3= v2;
691                                 n2->v1= v2;
692                         }
693                         else {
694                                 n1->v3= v1;
695                                 n2->v1= v1;
696                         }
697                 }
698                 else {
699                         n1->v3= n2->v1= mallocVert();
700                         n1->v3[0]= 0.5*(node->v1[0]+ node->v3[0]);
701                         n1->v3[1]= 0.5*(node->v1[1]+ node->v3[1]);
702                         n1->v3[2]= 0.5*(node->v1[2]+ node->v3[2]);
703                         n1->v3[3]= node->v3[3]; /* color */
704                 }
705         }
706         n1->area= AreaT3Dfl(n1->v1, n1->v2, n1->v3);
707         n2->area= AreaT3Dfl(n2->v1, n2->v2, n2->v3);
708
709 }
710
711
712 void subdivideNode(RNode *node, RNode *edge)
713 {
714         RNode *n1, *n2, *up;
715         float fu, fv, *v1, *v2;/*, AreaQ3Dfl(); ... from arithb... */
716         int uvl;
717         
718         if(Ntotnode>RG.maxnode) return;
719         
720         if(node->type==3) {
721                 subdivideTriNode(node, edge);
722                 return;
723         }
724
725         if(node->down1 || node->down2) {
726                 /* printf("subdivide Node: already done \n"); */
727                 return;
728         }
729         
730         /* defines subdivide direction */
731
732         if(edge==0) {
733                 /* areathreshold */
734                 if(node->area<nodelimit) {
735                         return;
736                 }
737                 fu= fabs(node->v1[0]- node->v2[0])+ fabs(node->v1[1]- node->v2[1]) +fabs(node->v1[2]- node->v2[2]);
738                 fv= fabs(node->v1[0]- node->v4[0])+ fabs(node->v1[1]- node->v4[1]) +fabs(node->v1[2]- node->v4[2]);
739                 if(fu>fv) uvl= 1;
740                 else uvl= 2;
741         }
742         else {
743                 if(edge==node->ed1 || edge==node->ed3) uvl= 1;
744                 else uvl= 2;
745         }
746         
747         /*  do neighbour nodes have to be deeper? Recursive! */
748         n1= n2= 0; 
749         if(uvl==1) {
750                 if(node->ed1 && node->ed1->down1==0) n1= node->ed1;
751                 if(node->ed3 && node->ed3->down1==0) n2= node->ed3;
752         }
753         else {
754                 if(node->ed2 && node->ed2->down1==0) n1= node->ed2;
755                 if(node->ed4 && node->ed4->down1==0) n2= node->ed4;
756         }
757         if(n1) {
758                 up= node->up;
759                 while(up) {
760                         if(n1->ed1==up || n1->ed2==up || n1->ed3==up || n1->ed4==up) {
761                                 /* printf("recurs subd\n"); */
762                                 subdivideNode(n1, up);
763                                 break;
764                         }
765                         up= up->up;
766                 }
767         }
768         if(n2) {
769                 up= node->up;
770                 while(up) {
771                         if(n2->ed1==up || n2->ed2==up || n2->ed3==up || n2->ed4==up) {
772                                 /* printf("recurs subd\n"); */
773                                 subdivideNode(n2, up);
774                                 break;
775                         }
776                         up= up->up;
777                 }
778         }
779
780         /* the subdividing */
781         n1= mallocNode();
782         memcpy(n1, node, sizeof(RNode));
783         n2= mallocNode();
784         memcpy(n2, node, sizeof(RNode));
785
786         n1->up= node;
787         n2->up= node;
788         
789         node->down1= n1;
790         node->down2= n2;
791
792         /* subdivide edge 1 and 3 */
793         if(uvl==1) {
794                 
795                 /* FIRST NODE  gets edge 2 */
796                 n1->ed4= n2;
797                 n1->lev4= 0;
798                 replaceAllNode(n1->ed2, n1);
799                 n1->lev1++;
800                 n1->lev3++;
801                 replaceTestNode(n1->ed1, &(n1->ed1), n1, n1->lev1, n1->v2);
802                 replaceTestNode(n1->ed3, &(n1->ed3), n1, n1->lev3, n1->v3);
803
804                 /* SECOND NODE gets edge 4 */
805                 n2->ed2= n1;
806                 n2->lev2= 0;
807                 replaceAllNode(n2->ed4, n2);
808                 n2->lev1++;
809                 n2->lev3++;
810                 replaceTestNode(n2->ed1, &(n2->ed1), n2, n2->lev1, n2->v1);
811                 replaceTestNode(n2->ed3, &(n2->ed3), n2, n2->lev3, n2->v4);
812                 
813                 /* NEW VERTEX from edge 1 */
814                 if( setvertexpointersNode(n1->ed1, n1, n1->lev1, &v1, &v2) ) {  /* nodes have equal levels */
815                         if(v1== n1->v2) {
816                                 n1->v1= v2;
817                                 n2->v2= v2;
818                         }
819                         else {
820                                 n1->v1= v1;
821                                 n2->v2= v1;
822                         }
823                 }
824                 else {
825                         n1->v1= n2->v2= mallocVert();
826                         n1->v1[0]= 0.5*(node->v1[0]+ node->v2[0]);
827                         n1->v1[1]= 0.5*(node->v1[1]+ node->v2[1]);
828                         n1->v1[2]= 0.5*(node->v1[2]+ node->v2[2]);
829                         n1->v1[3]= node->v1[3]; /* color */
830                 }
831                 
832                 /* NEW VERTEX from edge 3 */
833                 if( setvertexpointersNode(n1->ed3, n1, n1->lev3, &v1, &v2) ) {  /* nodes have equal levels */
834                         if(v1== n1->v3) {
835                                 n1->v4= v2;
836                                 n2->v3= v2;
837                         }
838                         else {
839                                 n1->v4= v1;
840                                 n2->v3= v1;
841                         }
842                 }
843                 else {
844                         n1->v4= n2->v3= mallocVert();
845                         n1->v4[0]= 0.5*(node->v3[0]+ node->v4[0]);
846                         n1->v4[1]= 0.5*(node->v3[1]+ node->v4[1]);
847                         n1->v4[2]= 0.5*(node->v3[2]+ node->v4[2]);
848                         n1->v4[3]= node->v4[3]; /* color */
849                 }
850         }
851         /* subdivide edge 2 and 4 */
852         else if(uvl==2) {
853                 
854                 /* FIRST NODE gets edge 1 */
855                 n1->ed3= n2;
856                 n1->lev3= 0;
857                 replaceAllNode(n1->ed1, n1);
858                 n1->lev2++;
859                 n1->lev4++;
860                 replaceTestNode(n1->ed2, &(n1->ed2), n1, n1->lev2, n1->v2);
861                 replaceTestNode(n1->ed4, &(n1->ed4), n1, n1->lev4, n1->v1);
862
863                 /* SECOND NODE gets edge 3 */
864                 n2->ed1= n1;
865                 n2->lev1= 0;
866                 replaceAllNode(n2->ed3, n2);
867                 n2->lev2++;
868                 n2->lev4++;
869                 replaceTestNode(n2->ed2, &(n2->ed2), n2, n2->lev2, n2->v3);
870                 replaceTestNode(n2->ed4, &(n2->ed4), n2, n2->lev4, n2->v4);
871
872                 /* NEW VERTEX from edge 2 */
873                 if( setvertexpointersNode(n1->ed2, n1, n1->lev2, &v1, &v2) ) {  /* nodes have equal levels */
874                         if(v1== n1->v2) {
875                                 n1->v3= v2;
876                                 n2->v2= v2;
877                         }
878                         else {
879                                 n1->v3= v1;
880                                 n2->v2= v1;
881                         }
882                 }
883                 else {
884                         n1->v3= n2->v2= mallocVert();
885                         n1->v3[0]= 0.5*(node->v2[0]+ node->v3[0]);
886                         n1->v3[1]= 0.5*(node->v2[1]+ node->v3[1]);
887                         n1->v3[2]= 0.5*(node->v2[2]+ node->v3[2]);
888                         n1->v3[3]= node->v3[3]; /* color */
889                 }
890
891                 /* NEW VERTEX from edge 4 */
892                 if( setvertexpointersNode(n1->ed4, n1, n1->lev4, &v1, &v2) ) {  /* nodes have equal levels */
893                         if(v1== n1->v1) {
894                                 n1->v4= v2;
895                                 n2->v1= v2;
896                         }
897                         else {
898                                 n1->v4= v1;
899                                 n2->v1= v1;
900                         }
901                 }
902                 else {
903                         n1->v4= n2->v1= mallocVert();
904                         n1->v4[0]= 0.5*(node->v1[0]+ node->v4[0]);
905                         n1->v4[1]= 0.5*(node->v1[1]+ node->v4[1]);
906                         n1->v4[2]= 0.5*(node->v1[2]+ node->v4[2]);
907                         n1->v4[3]= node->v4[3]; /* color */
908                 }
909         }
910         
911         n1->area= AreaQ3Dfl(n1->v1, n1->v2, n1->v3, n1->v4);
912         n2->area= AreaQ3Dfl(n2->v1, n2->v2, n2->v3, n2->v4);
913
914 }
915
916 int comparelevel(RNode *node, RNode *nb, int level)
917 {
918         /* recursive descent: test with deepest node */
919         /* return 1 means equal or higher */
920         
921         if(nb==0) return 1;
922         
923         if(nb->down1) {
924                 return 0;
925                 
926                 /*              THERE IS AN ERROR HERE, BUT WHAT?  (without this function the system
927                         works too, but is slower) (ton) */
928
929                 /*
930                 n1= nb->down1;
931                 if(n1->ed1==node) return comparelevel(node, n1, level);
932                 if(n1->ed2==node) return comparelevel(node, n1, level);
933                 if(n1->ed3==node) return comparelevel(node, n1, level);
934                 if(n1->ed4==node) return comparelevel(node, n1, level);
935                 n1= nb->down2;
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                 printf(" dit kan niet ");
941                 return 0;
942                 */
943                 
944         }
945         
946         if(nb->down1==0) {
947                 /* if(nb->ed1==node) return (nb->lev1<=level); */
948                 /* if(nb->ed2==node) return (nb->lev2<=level); */
949                 /* if(nb->ed3==node) return (nb->lev3<=level); */
950                 /* if(nb->ed4==node) return (nb->lev4<=level); */
951                 
952                 return 1;       /* is higher node */
953         }
954         return 1;
955 }
956
957 void deleteTriNodes(RNode *node)        /* both children of node */
958 {
959         RNode *n1, *n2;
960         
961         /* if neighbour nodes are deeper: no delete */
962         /* just test 2 nodes, from the others the level doesn't change */
963         
964         n1= node->down1;
965         n2= node->down2;
966
967         if(n1==0 || n2==0) return;
968         
969         if(n1->down1 || n2->down1) return;
970         
971         /* at the edges no subdivided node is allowed */
972
973         if(n1->ed1 && n1->ed1->down1) return;
974         if(n1->ed2 && n1->ed2->down1) return;
975         if(n1->ed3 && n1->ed3->down1) return;
976
977         if(n2->ed1 && n2->ed1->down1) return;
978         if(n2->ed2 && n2->ed2->down1) return;
979         if(n2->ed3 && n2->ed3->down1) return;
980                                 
981         replaceAllNodeInv(n1->ed1, n1);
982         replaceAllNodeInv(n1->ed2, n1);
983         replaceAllNodeInv(n1->ed3, n1);
984
985         replaceAllNodeUp(n1->ed1, n1);
986         replaceAllNodeUp(n1->ed2, n1);
987         replaceAllNodeUp(n1->ed3, n1);
988
989         replaceAllNodeInv(n2->ed1, n2);
990         replaceAllNodeInv(n2->ed2, n2);
991         replaceAllNodeInv(n2->ed3, n2);
992
993         replaceAllNodeUp(n2->ed1, n2);
994         replaceAllNodeUp(n2->ed2, n2);
995         replaceAllNodeUp(n2->ed3, n2);
996
997         n1->down1= (RNode *)12; /* for debug */
998         n2->down1= (RNode *)12;
999         
1000         freeNode(n1);
1001         freeNode(n2);
1002         node->down1= node->down2= 0;
1003         
1004 }
1005
1006         /* both children of node */
1007 void deleteNodes(RNode *node)
1008 {
1009         RNode *n1, *n2;
1010
1011         /* if neighbour nodes are deeper: no delete */
1012         /* just test 2 nodes, from the others the level doesn't change */
1013
1014         if(node->type==3) {
1015                 deleteTriNodes(node);
1016                 return;
1017         }
1018         
1019         n1= node->down1;
1020         n2= node->down2;
1021
1022         if(n1==0 || n2==0) return;
1023         
1024         if(n1->down1 || n2->down1) return;
1025         
1026         if(n1->ed3==n2) {
1027
1028                 /* at the edges no subdivided node is allowed */
1029
1030                 if(n1->ed1 && n1->ed1->down1) return;
1031                 if(n1->ed2 && n1->ed2->down1) return;
1032                 if(n1->ed4 && n1->ed4->down1) return;
1033
1034                 if(n2->ed2 && n2->ed2->down1) return;
1035                 if(n2->ed3 && n2->ed3->down1) return;
1036                 if(n2->ed4 && n2->ed4->down1) return;
1037                                         
1038                 replaceAllNodeInv(n1->ed1, n1);
1039                 replaceAllNodeInv(n1->ed2, n1);
1040                 replaceAllNodeInv(n1->ed4, n1);
1041
1042                 replaceAllNodeUp(n1->ed1, n1);
1043                 replaceAllNodeUp(n1->ed2, n1);
1044                 replaceAllNodeUp(n1->ed4, n1);
1045
1046                 replaceAllNodeInv(n2->ed2, n2);
1047                 replaceAllNodeInv(n2->ed3, n2);
1048                 replaceAllNodeInv(n2->ed4, n2);
1049
1050                 replaceAllNodeUp(n2->ed2, n2);
1051                 replaceAllNodeUp(n2->ed3, n2);
1052                 replaceAllNodeUp(n2->ed4, n2);
1053
1054                 n1->down1= (RNode *)12; /* for debug */
1055                 n2->down1= (RNode *)12;
1056                 
1057                 freeNode(n1);
1058                 freeNode(n2);
1059                 node->down1= node->down2= 0;
1060                 
1061                 return;
1062         }
1063         else if(n1->ed4==n2) {
1064
1065                 if(n1->ed1 && n1->ed1->down1) return;
1066                 if(n1->ed2 && n1->ed2->down1) return;
1067                 if(n1->ed3 && n1->ed3->down1) return;
1068
1069                 if(n2->ed1 && n2->ed1->down1) return;
1070                 if(n2->ed3 && n2->ed3->down1) return;
1071                 if(n2->ed4 && n2->ed4->down1) return;
1072                                         
1073                 replaceAllNodeInv(n1->ed1, n1);
1074                 replaceAllNodeInv(n1->ed2, n1);
1075                 replaceAllNodeInv(n1->ed3, n1);
1076
1077                 replaceAllNodeUp(n1->ed1, n1);
1078                 replaceAllNodeUp(n1->ed2, n1);
1079                 replaceAllNodeUp(n1->ed3, n1);
1080
1081                 replaceAllNodeInv(n2->ed1, n2);
1082                 replaceAllNodeInv(n2->ed3, n2);
1083                 replaceAllNodeInv(n2->ed4, n2);
1084
1085                 replaceAllNodeUp(n2->ed1, n2);
1086                 replaceAllNodeUp(n2->ed3, n2);
1087                 replaceAllNodeUp(n2->ed4, n2);
1088
1089                 n1->down1= (RNode *)12; /* for debug */
1090                 n2->down1= (RNode *)12;
1091                 
1092                 freeNode(n1);
1093                 freeNode(n2);
1094                 node->down1= node->down2= 0;
1095                 
1096                 return;
1097         }
1098
1099 }
1100
1101