updated .c files to include:
[blender.git] / source / blender / radiosity / intern / source / radfactors.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     formfactors.c       nov/dec 1992
34
35     $Id$
36
37  *************************************** */
38
39 #include <stdlib.h>
40 #include <string.h>
41 #include <math.h>
42
43 #include "MEM_guardedalloc.h"
44
45 #include "BLI_blenlib.h"
46 #include "BLI_arithb.h"
47 #include "BLI_rand.h"
48
49 #include "BKE_utildefines.h"
50 #include "BKE_global.h"
51 #include "BKE_main.h"
52
53 #include "BIF_screen.h"
54
55 #include "radio.h"
56 #include "render.h"       /* for `RE_zbufferall_radio and RE_zbufferall_radio */
57 #include "blendertimer.h" /* timer functions */
58
59 #ifdef HAVE_CONFIG_H
60 #include <config.h>
61 #endif
62
63 /* locals */
64 void rad_setmatrices(RadView *vw);
65 void clearsubflagelem(RNode *rn);
66 void setsubflagelem(RNode *rn);
67
68 RadView hemitop, hemiside;
69
70 float calcStokefactor(RPatch *shoot, RPatch *rp, RNode *rn, float *area)
71 {
72         float tvec[3], fac;
73         float vec[4][3];        /* vectoren van shootcent naar vertices rp */
74         float cross[4][3];      /* uitprodukten hiervan */
75         float rad[4];   /* hoeken tussen vecs */
76
77         /* test op richting */
78         VecSubf(tvec, shoot->cent, rp->cent);
79         if( tvec[0]*shoot->norm[0]+ tvec[1]*shoot->norm[1]+ tvec[2]*shoot->norm[2]>0.0)
80                 return 0.0;
81         
82         if(rp->type==4) {
83
84                 /* hoekvectors */
85                 VecSubf(vec[0], shoot->cent, rn->v1);
86                 VecSubf(vec[1], shoot->cent, rn->v2);
87                 VecSubf(vec[2], shoot->cent, rn->v3);
88                 VecSubf(vec[3], shoot->cent, rn->v4);
89
90                 Normalise(vec[0]);
91                 Normalise(vec[1]);
92                 Normalise(vec[2]);
93                 Normalise(vec[3]);
94
95                 /* uitprod */
96                 Crossf(cross[0], vec[0], vec[1]);
97                 Crossf(cross[1], vec[1], vec[2]);
98                 Crossf(cross[2], vec[2], vec[3]);
99                 Crossf(cross[3], vec[3], vec[0]);
100                 Normalise(cross[0]);
101                 Normalise(cross[1]);
102                 Normalise(cross[2]);
103                 Normalise(cross[3]);
104
105                 /* hoeken */
106                 rad[0]= vec[0][0]*vec[1][0]+ vec[0][1]*vec[1][1]+ vec[0][2]*vec[1][2];
107                 rad[1]= vec[1][0]*vec[2][0]+ vec[1][1]*vec[2][1]+ vec[1][2]*vec[2][2];
108                 rad[2]= vec[2][0]*vec[3][0]+ vec[2][1]*vec[3][1]+ vec[2][2]*vec[3][2];
109                 rad[3]= vec[3][0]*vec[0][0]+ vec[3][1]*vec[0][1]+ vec[3][2]*vec[0][2];
110
111                 rad[0]= acos(rad[0]);
112                 rad[1]= acos(rad[1]);
113                 rad[2]= acos(rad[2]);
114                 rad[3]= acos(rad[3]);
115
116                 /* Stoke formule */
117                 VecMulf(cross[0], rad[0]);
118                 VecMulf(cross[1], rad[1]);
119                 VecMulf(cross[2], rad[2]);
120                 VecMulf(cross[3], rad[3]);
121
122                 VECCOPY(tvec, shoot->norm);
123                 fac=  tvec[0]*cross[0][0]+ tvec[1]*cross[0][1]+ tvec[2]*cross[0][2];
124                 fac+= tvec[0]*cross[1][0]+ tvec[1]*cross[1][1]+ tvec[2]*cross[1][2];
125                 fac+= tvec[0]*cross[2][0]+ tvec[1]*cross[2][1]+ tvec[2]*cross[2][2];
126                 fac+= tvec[0]*cross[3][0]+ tvec[1]*cross[3][1]+ tvec[2]*cross[3][2];
127         }
128         else {
129                 /* hoekvectors */
130                 VecSubf(vec[0], shoot->cent, rn->v1);
131                 VecSubf(vec[1], shoot->cent, rn->v2);
132                 VecSubf(vec[2], shoot->cent, rn->v3);
133
134                 Normalise(vec[0]);
135                 Normalise(vec[1]);
136                 Normalise(vec[2]);
137
138                 /* uitprod */
139                 Crossf(cross[0], vec[0], vec[1]);
140                 Crossf(cross[1], vec[1], vec[2]);
141                 Crossf(cross[2], vec[2], vec[0]);
142                 Normalise(cross[0]);
143                 Normalise(cross[1]);
144                 Normalise(cross[2]);
145
146                 /* hoeken */
147                 rad[0]= vec[0][0]*vec[1][0]+ vec[0][1]*vec[1][1]+ vec[0][2]*vec[1][2];
148                 rad[1]= vec[1][0]*vec[2][0]+ vec[1][1]*vec[2][1]+ vec[1][2]*vec[2][2];
149                 rad[2]= vec[2][0]*vec[0][0]+ vec[2][1]*vec[0][1]+ vec[2][2]*vec[0][2];
150
151                 rad[0]= acos(rad[0]);
152                 rad[1]= acos(rad[1]);
153                 rad[2]= acos(rad[2]);
154
155                 /* Stoke formule */
156                 VecMulf(cross[0], rad[0]);
157                 VecMulf(cross[1], rad[1]);
158                 VecMulf(cross[2], rad[2]);
159
160                 VECCOPY(tvec, shoot->norm);
161                 fac=  tvec[0]*cross[0][0]+ tvec[1]*cross[0][1]+ tvec[2]*cross[0][2];
162                 fac+= tvec[0]*cross[1][0]+ tvec[1]*cross[1][1]+ tvec[2]*cross[1][2];
163                 fac+= tvec[0]*cross[2][0]+ tvec[1]*cross[2][1]+ tvec[2]*cross[2][2];
164         }
165
166         *area= -fac/(2.0*PI);
167         return (*area * (shoot->area/rn->area));
168 }
169
170
171 void calcTopfactors()
172 {
173         float xsq , ysq, xysq;
174         float n;
175         float *fp;
176         int a, b, hres;
177
178         fp = RG.topfactors;
179         hres= RG.hemires/2;
180         n= hres;
181
182         for (a=0; a<hres; a++) {
183
184                 ysq= (n- ((float)a+0.5))/n;
185                 ysq*= ysq;
186
187                 for ( b=0 ; b<hres ; b++ ) {
188
189                         xsq= ( n-((float)b+ 0.5) )/n;
190                         xsq*= xsq;
191                         xysq=  xsq+ ysq+ 1.0 ;
192                         xysq*= xysq;
193
194                         *fp++ = 1.0/(xysq* PI* n*n);
195                 }
196         }
197
198 }
199
200 void calcSidefactors()
201 {
202         float xsq , ysq, xysq;
203         float n, y;
204         float *fp;
205         int a, b, hres;
206
207         fp = RG.sidefactors;
208         hres= RG.hemires/2;
209         n= hres;
210
211         for (a=0; a<hres; a++) {
212
213                 y= (n- ((float)a+0.5))/n;
214                 ysq= y*y;
215
216                 for ( b=0 ; b<hres ; b++ ) {
217
218                         xsq= ( n-((float)b+ 0.5) )/n;
219                         xsq*= xsq;
220                         xysq=  xsq+ ysq+ 1.0 ;
221                         xysq*= xysq;
222
223                         *fp++ = y/(xysq* PI* n*n);
224                 }
225         }
226
227 }
228
229
230 void initradiosity()
231 {
232         /* alloceert en maakt LUTs voor top/side factors */
233         /* alloceert en maakt index array */
234         int a, hres2;
235
236         if(RG.topfactors) MEM_freeN(RG.topfactors);
237         if(RG.sidefactors) MEM_freeN(RG.sidefactors);
238         if(RG.index) MEM_freeN(RG.index);
239
240         RG.topfactors= MEM_callocN(RG.hemires*RG.hemires, "initradiosity");
241         calcTopfactors();
242         RG.sidefactors= MEM_callocN(RG.hemires*RG.hemires, "initradiosity1");
243         calcSidefactors();
244
245         RG.index= MEM_callocN(4*RG.hemires, "initradiosity3");
246         hres2= RG.hemires/2;
247         for(a=0; a<RG.hemires; a++) {
248                 RG.index[a]=  a<hres2 ? a: (hres2-1-( a % hres2 ));
249         }
250
251 }
252
253 void rad_make_hocos(RadView *vw)
254 {
255         /* float hoco[4]; */
256         /* int a; */
257         
258         /* for(a=0; a< R.totvert;a++) { */
259         /*      projectvert(vec, ver->ho); */
260         /*      ver->clip = testclip(ver->ho); */
261 /*  */
262         /* } */
263 }
264
265 void rad_setmatrices(RadView *vw)           /* voor hemi's */
266 {
267         float up1[3], len, twist;
268
269         i_lookat(vw->cam[0], vw->cam[1], vw->cam[2], vw->tar[0], vw->tar[1], vw->tar[2], 0, vw->viewmat);
270         up1[0] = vw->viewmat[0][0]*vw->up[0] + vw->viewmat[1][0]*vw->up[1] + vw->viewmat[2][0]*vw->up[2];
271         up1[1] = vw->viewmat[0][1]*vw->up[0] + vw->viewmat[1][1]*vw->up[1] + vw->viewmat[2][1]*vw->up[2];
272         up1[2] = vw->viewmat[0][2]*vw->up[0] + vw->viewmat[1][2]*vw->up[1] + vw->viewmat[2][2]*vw->up[2];
273
274         len= up1[0]*up1[0]+up1[1]*up1[1];
275         if(len>0.0) {
276                 twist= -atan2(up1[0], up1[1]);
277         }
278         else twist= 0.0;
279         
280         i_lookat(vw->cam[0], vw->cam[1], vw->cam[2], vw->tar[0], vw->tar[1], vw->tar[2], (180.0*twist/M_PI), vw->viewmat);
281
282         /* window matrix was set in inithemiwindows */
283         
284 }
285
286
287 void hemizbuf(RadView *vw)
288 {
289         float *factors;
290         unsigned int *rz;
291         int a, b, inda, hres;
292
293         rad_setmatrices(vw);
294         RE_zbufferall_radio(vw, RG.elem, RG.totelem);
295
296         /* factors tellen */
297         if(vw->recty==vw->rectx) factors= RG.topfactors;
298         else factors= RG.sidefactors;
299         hres= RG.hemires/2;
300
301         rz= vw->rect;
302         for(a=0; a<vw->recty; a++) {
303                 inda= hres*RG.index[a];
304                 for(b=0; b<vw->rectx; b++, rz++) {
305                         if(*rz<RG.totelem) {
306                                 RG.formfactors[*rz]+= factors[inda+RG.index[b]];
307                         }
308                 }
309         }
310 }
311
312 int makeformfactors(RPatch *shoot)
313 {
314         RNode **re;
315         float len, vec[3], up[3], side[3], tar[5][3], *fp;
316         int a, overfl;
317
318         if(RG.totelem==0) return 0;
319
320         memset(RG.formfactors, 0, 4*RG.totelem);
321
322         /* set up hemiview */
323         /* first: random upvector */
324         do {
325                 vec[0]= (float)BLI_drand();
326                 vec[1]= (float)BLI_drand();
327                 vec[2]= (float)BLI_drand();
328                 Crossf(up, shoot->norm, vec);
329                 len= Normalise(up);
330         } while(len==0.0 || len>1.0);
331
332         VECCOPY(hemitop.up, up);
333         VECCOPY(hemiside.up, shoot->norm);
334
335         Crossf(side, shoot->norm, up);
336
337         /* five targets */
338         VecAddf(tar[0], shoot->cent, shoot->norm);
339         VecAddf(tar[1], shoot->cent, up);
340         VecSubf(tar[2], shoot->cent, up);
341         VecAddf(tar[3], shoot->cent, side);
342         VecSubf(tar[4], shoot->cent, side);
343
344         /* camera */
345         VECCOPY(hemiside.cam, shoot->cent);
346         VECCOPY(hemitop.cam, shoot->cent);
347
348         /* do it! */
349         VECCOPY(hemitop.tar, tar[0]);
350         hemizbuf(&hemitop);
351
352         for(a=1; a<5; a++) {
353                 VECCOPY(hemiside.tar, tar[a]);
354                 hemizbuf(&hemiside);
355         }
356
357         /* convert factors to real radiosity */
358         re= RG.elem;
359         fp= RG.formfactors;
360
361         overfl= 0;
362         for(a= RG.totelem; a>0; a--, re++, fp++) {
363         
364                 if(*fp!=0.0) {
365                         
366                         *fp *= shoot->area/(*re)->area;
367
368                         if(*fp>1.0) {
369                                 overfl= 1;
370                                 *fp= 1.0001;
371                         }
372                 }
373         }
374         
375         if(overfl) {
376                 /*
377                 drawOverflowElem();
378                 while(get_mbut()&L_MOUSE==0) {
379                         if(get_mbut()&M_MOUSE) {
380                                 viewmove();
381                                 drawpatch_ext(shoot,0xFF77FF);
382                                 drawOverflowElem();
383                         }
384                 }
385                 */
386                 if(shoot->first->down1) {
387                         splitpatch(shoot);
388                         return 0;
389                 }
390         }
391         
392         return 1;
393 }
394
395 void applyformfactors(RPatch *shoot)
396 {
397         RPatch *rp;
398         RNode **el, *rn;
399         float *fp, *ref, unr, ung, unb, r, g, b, w;
400         int a;
401
402         unr= shoot->unshot[0];
403         ung= shoot->unshot[1];
404         unb= shoot->unshot[2];
405
406         fp= RG.formfactors;
407         el= RG.elem;
408         for(a=0; a<RG.totelem; a++, el++, fp++) {
409                 rn= *el;
410                 if(*fp!= 0.0) {
411                         rp= rn->par;
412                         ref= rp->ref;
413                         
414                         r= (*fp)*unr*ref[0];
415                         g= (*fp)*ung*ref[1];
416                         b= (*fp)*unb*ref[2];
417
418                         w= rn->area/rp->area;
419                         rn->totrad[0]+= r;
420                         rn->totrad[1]+= g;
421                         rn->totrad[2]+= b;
422                         
423                         rp->unshot[0]+= w*r;
424                         rp->unshot[1]+= w*g;
425                         rp->unshot[2]+= w*b;
426                 }
427         }
428         
429         shoot->unshot[0]= shoot->unshot[1]= shoot->unshot[2]= 0.0;
430 }
431
432 RPatch *findshootpatch()
433 {
434         RPatch *rp, *shoot;
435         float energy, maxenergy;
436
437         shoot= 0;
438         maxenergy= 0.0;
439         rp= RG.patchbase.first;
440         while(rp) {
441                 energy= rp->unshot[0]*rp->area;
442                 energy+= rp->unshot[1]*rp->area;
443                 energy+= rp->unshot[2]*rp->area;
444
445                 if(energy>maxenergy) {
446                         shoot= rp;
447                         maxenergy= energy;
448                 }
449                 rp= rp->next;
450         }
451
452         if(shoot) {
453                 maxenergy/= RG.totenergy;
454                 if(maxenergy<RG.convergence) return 0;
455         }
456
457         return shoot;
458 }
459
460 void setnodeflags(RNode *rn, int flag, int set)
461 {
462         
463         if(rn->down1) {
464                 setnodeflags(rn->down1, flag, set);
465                 setnodeflags(rn->down2, flag, set);
466         }
467         else {
468                 if(set) rn->f |= flag;
469                 else rn->f &= ~flag;
470         }
471 }
472
473 void backface_test(RPatch *shoot)
474 {
475         RPatch *rp;
476         float tvec[3];
477         
478         rp= RG.patchbase.first;
479         while(rp) {
480                 if(rp!=shoot) {
481                 
482                         VecSubf(tvec, shoot->cent, rp->cent);
483                         if( tvec[0]*shoot->norm[0]+ tvec[1]*shoot->norm[1]+ tvec[2]*shoot->norm[2]>0.0) {
484                                 setnodeflags(rp->first, RAD_BACKFACE, 1);
485                         }
486                 }
487                 rp= rp->next;
488         }
489 }
490
491 void clear_backface_test()
492 {
493         RNode **re;
494         int a;
495         
496         re= RG.elem;
497         for(a= RG.totelem-1; a>=0; a--, re++) {
498                 (*re)->f &= ~RAD_BACKFACE;
499         }
500         
501 }
502
503 void rad_init_energy()
504 {
505         /* call before shooting */
506         /* keep patches and elements, clear all data */
507         RNode **el, *rn;
508         RPatch *rp;
509         int a;
510
511         el= RG.elem;
512         for(a=RG.totelem; a>0; a--, el++) {
513                 rn= *el;
514                 VECCOPY(rn->totrad, rn->par->emit);
515         }
516
517         RG.totenergy= 0.0;
518         rp= RG.patchbase.first;
519         while(rp) {
520                 VECCOPY(rp->unshot, rp->emit);
521
522                 RG.totenergy+= rp->unshot[0]*rp->area;
523                 RG.totenergy+= rp->unshot[1]*rp->area;
524                 RG.totenergy+= rp->unshot[2]*rp->area;
525
526                 rp->f= 0;
527                 
528                 rp= rp->next;
529         }
530 }
531
532 void progressiverad()
533 {
534         RPatch *shoot;
535         int it= 0;
536
537         rad_printstatus();
538         rad_init_energy();
539
540         shoot=findshootpatch();
541
542         while( shoot ) {
543
544                 setnodeflags(shoot->first, RAD_SHOOT, 1);
545                 
546                 backface_test(shoot);
547                 
548                 drawpatch_ext(shoot, 0x88FF00);
549
550                 if( makeformfactors(shoot) ) {
551                 
552                         applyformfactors(shoot);
553         
554                         it++;
555                         set_timecursor(it);
556                         if( (it & 3)==1 ) {
557                                 make_node_display();
558                                 rad_forcedraw();
559                         }
560                         setnodeflags(shoot->first, RAD_SHOOT, 0);
561                 }
562                 
563                 clear_backface_test();
564                 
565                 if(MISC_test_break()) break;
566                 if(RG.maxiter && RG.maxiter<=it) break;
567
568                 shoot=findshootpatch();
569
570         }
571         
572 }
573
574
575 /* *************  subdivideshoot *********** */
576
577 void minmaxradelem(RNode *rn, float *min, float *max)
578 {
579         int c;
580         
581         if(rn->down1) {
582                 minmaxradelem(rn->down1, min, max);
583                 minmaxradelem(rn->down2, min, max);
584         }
585         else {
586                 for(c=0; c<3; c++) {
587                         min[c]= MIN2(min[c], rn->totrad[c]);
588                         max[c]= MAX2(max[c], rn->totrad[c]);
589                 }
590         }
591 }
592
593 void minmaxradelemfilt(RNode *rn, float *min, float *max, float *errmin, float *errmax)
594 {
595         float col[3], area;
596         int c;
597         
598         if(rn->down1) {
599                 minmaxradelemfilt(rn->down1, min, max, errmin, errmax);
600                 minmaxradelemfilt(rn->down2, min, max, errmin, errmax);
601         }
602         else {
603                 VECCOPY(col, rn->totrad);
604
605                 for(c=0; c<3; c++) {
606                         min[c]= MIN2(min[c], col[c]);
607                         max[c]= MAX2(max[c], col[c]);
608                 }
609
610                 VecMulf(col, 2.0);
611                 area= 2.0;
612                 if(rn->ed1) {
613                         VecAddf(col, rn->ed1->totrad, col);
614                         area+= 1.0;
615                 }
616                 if(rn->ed2) {
617                         VecAddf(col, rn->ed2->totrad, col);
618                         area+= 1.0;
619                 }
620                 if(rn->ed3) {
621                         VecAddf(col, rn->ed3->totrad, col);
622                         area+= 1.0;
623                 }
624                 if(rn->ed4) {
625                         VecAddf(col, rn->ed4->totrad, col);
626                         area+= 1.0;
627                 }
628                 VecMulf(col, 1.0/area);
629
630                 for(c=0; c<3; c++) {
631                         errmin[c]= MIN2(errmin[c], col[c]);
632                         errmax[c]= MAX2(errmax[c], col[c]);
633                 }
634         }
635 }
636
637 void setsubflagelem(RNode *rn)
638 {
639         
640         if(rn->down1) {
641                 setsubflagelem(rn->down1);
642                 setsubflagelem(rn->down2);
643         }
644         else {
645                 rn->f |= RAD_SUBDIV;
646         }
647 }
648
649 void clearsubflagelem(RNode *rn)
650 {
651         
652         if(rn->down1) {
653                 setsubflagelem(rn->down1);
654                 setsubflagelem(rn->down2);
655         }
656         else {
657                 rn->f &= ~RAD_SUBDIV;
658         }
659 }
660
661 void subdivideshootElements(int it)
662 {
663         RPatch *rp, *shoot;
664         RNode **el, *rn;
665         float *fp, err, stoke, area, min[3], max[3], errmin[3], errmax[3];
666         int a, b, c, d, e, f, contin;
667         int maxlamp;
668         
669         if(RG.maxsublamp==0) maxlamp= RG.totlamp;
670         else maxlamp= RG.maxsublamp;
671         
672         while(it) {
673                 rad_printstatus();
674                 rad_init_energy();
675                 it--;
676                 
677                 for(a=0; a<maxlamp; a++) {
678                         shoot= findshootpatch();
679                         if(shoot==0) break;
680                         
681                         set_timecursor(a);
682                         drawpatch_ext(shoot, 0x88FF00);
683                         
684                         setnodeflags(shoot->first, RAD_SHOOT, 1);
685                         if( makeformfactors(shoot) ) {
686                         
687                                 fp= RG.formfactors;
688                                 el= RG.elem;
689                                 for(b=RG.totelem; b>0; b--, el++) {
690                                         rn= *el;
691                                         
692                                         if( (rn->f & RAD_SUBDIV)==0 && *fp!=0.0) {
693                                                 if(rn->par->emit[0]+rn->par->emit[1]+rn->par->emit[2]==0.0) {
694                                                 
695                                                         stoke= calcStokefactor(shoot, rn->par, rn, &area);
696                                                         if(stoke!= 0.0) {
697                                                         
698                                                                 err= *fp/stoke;
699                                                                 
700                                                                 /* area error */
701                                                                 area*=(0.5*RG.hemires*RG.hemires);
702                                                                 
703                                                                 if(area>35.0) {
704                                                                         if(err<0.95 || err>1.05) {
705                                                                                 if(err>0.05) {
706                                                                                         rn->f |= RAD_SUBDIV;
707                                                                                         rn->par->f |= RAD_SUBDIV;
708                                                                                 }
709                                                                         }
710                                                                 }
711                                                         }
712                                                         
713                                                 }
714                                         }
715                                         
716                                         fp++;
717         
718                                 }
719         
720                                 applyformfactors(shoot);
721                 
722                                 if( (a & 3)==1 ) {
723                                         make_node_display();
724                                         rad_forcedraw();
725                                 }
726                                 
727                                 setnodeflags(shoot->first, RAD_SHOOT, 0);
728                         }
729                         else a--;
730                         
731                         if(MISC_test_break()) break;
732                 }
733                 
734                 /* test op extreem weinig kleurverloop binnen patch met subdivflag */
735                 
736                 rp= RG.patchbase.first;
737                 
738                 while(rp) {
739                         if(rp->f & RAD_SUBDIV) {                /* rp heeft elems die moet gesubd */
740                                 /* minstens 4 levels diep */
741                                 rn= rp->first->down1;
742                                 if(rn) {
743                                         rn= rn->down1;
744                                         if(rn) {
745                                                 rn= rn->down1;
746                                                 if(rn) rn= rn->down1;
747                                         }
748                                 }
749                                 if(rn) {
750                                         min[0]= min[1]= min[2]= 1.0e10;
751                                         max[0]= max[1]= max[2]= -1.0e10;
752                                         /* errmin en max zijn de gefilterde kleuren */
753                                         errmin[0]= errmin[1]= errmin[2]= 1.0e10;
754                                         errmax[0]= errmax[1]= errmax[2]= -1.0e10;
755                                         minmaxradelemfilt(rp->first, min, max, errmin, errmax);
756                                         
757                                         /* verschil tussen kleuren klein: toch maar niet subd */
758                                         /* test ook voor de gefilterde: maar kritischer */
759                                         
760                                         contin= 0;
761                                         a= abs( calculatecolor(min[0])-calculatecolor(max[0]));
762                                         b= abs( calculatecolor(errmin[0])-calculatecolor(errmax[0]));
763                                         if(a<15 || b<7) {
764                                                 c= abs( calculatecolor(min[1])-calculatecolor(max[1]));
765                                                 d= abs( calculatecolor(errmin[1])-calculatecolor(errmax[1]));
766                                                 if(c<15 || d<7) {
767                                                         e= abs( calculatecolor(min[2])-calculatecolor(max[2]));
768                                                         f= abs( calculatecolor(errmin[2])-calculatecolor(errmax[2]));
769                                                         if(e<15 || f<7) {
770                                                                 contin= 1;
771                                                                 clearsubflagelem(rp->first);
772                                                                 /* printf("%d %d %d %d %d %d\n", a, b, c, d, e, f); */
773                                                         }
774                                                 }
775                                         }
776                                         if(contin) {
777                                                 drawpatch_ext(rp, 0xFFFF);
778                                         }
779                                 }
780                         }
781                         rp->f &= ~RAD_SUBDIV;
782                         rp= rp->next;
783                 }
784                 
785                 contin= 0;
786                 
787                 el= RG.elem;
788                 for(b=RG.totelem; b>0; b--, el++) {
789                         rn= *el;
790                         if(rn->f & RAD_SUBDIV) {
791                                 rn->f-= RAD_SUBDIV;
792                                 subdivideNode(rn, 0);
793                                 if(rn->down1) {
794                                         subdivideNode(rn->down1, 0);
795                                         subdivideNode(rn->down2, 0);
796                                         contin= 1;
797                                 }
798                         }
799                 }
800                 makeGlobalElemArray();
801
802                 if(contin==0 || MISC_test_break()) break;
803         }
804         
805         make_node_display();
806 }
807
808 void subdivideshootPatches(int it)
809 {
810         RPatch *rp, *shoot, *next;
811         float *fp, err, stoke, area;
812         int a, contin;
813         int maxlamp;
814         
815         if(RG.maxsublamp==0) maxlamp= RG.totlamp;
816         else maxlamp= RG.maxsublamp;
817         
818         while(it) {
819                 rad_printstatus();
820                 rad_init_energy();
821                 it--;
822                 
823                 for(a=0; a<maxlamp; a++) {
824                         shoot= findshootpatch();
825                         if(shoot==0) break;
826                         
827                         set_timecursor(a);
828                         drawpatch_ext(shoot, 0x88FF00);
829                         
830                         setnodeflags(shoot->first, RAD_SHOOT, 1);
831                         
832                         if( makeformfactors(shoot) ) {
833                         
834                                 fp= RG.formfactors;
835                                 rp= RG.patchbase.first;
836                                 while(rp) {
837                                         if(*fp!=0.0 && rp!=shoot) {
838                                         
839                                                 stoke= calcStokefactor(shoot, rp, rp->first, &area);
840                                                 if(stoke!= 0.0) {
841                                                         if(area>.1) {   /* ontvangt patch meer dan (ong)10% van energie? */
842                                                                 rp->f= RAD_SUBDIV;
843                                                         }
844                                                         else {
845                                         
846                                                                 err= *fp/stoke;
847                                                         
848                                                                 /* area error */
849                                                                 area*=(0.5*RG.hemires*RG.hemires);
850                                                                 
851                                                                 if(area>45.0) {
852                                                                         if(err<0.95 || err>1.05) {
853                                                                                 if(err>0.05) {
854                                                                                                                 
855                                                                                         rp->f= RAD_SUBDIV;
856                                                                                         
857                                                                                         /* if(get_qual()&LR_SHIFTKEY);
858                                                                                         else {
859                                                                                                 drawpatch_ext(rp, 0xFF77FF);
860                                         
861                                                                                                 printf("Pa hemi %f stoke %f err %f area %f\n", *fp, stoke, err, area);
862                                         
863                                                                                                 while(get_mbut()&L_MOUSE==0);
864                                                                                                 while(get_mbut()&L_MOUSE);
865                                                                                         }
866                                                                                         */
867                                                                                 }
868                                                                         }
869                                                                 }
870                                                         }
871                                                 }
872                                         }
873                                         fp++;
874         
875                                         rp= rp->next;
876                                 }
877         
878                                 applyformfactors(shoot);
879                 
880                                 if( (a & 3)==1 ) {
881                                         make_node_display();
882                                         rad_forcedraw();
883                                 }
884                                 
885                                 setnodeflags(shoot->first, RAD_SHOOT, 0);
886                                 
887                                 if(MISC_test_break()) break;
888                         }
889                         else a--;
890                         
891                 }
892                 
893                 contin= 0;
894
895                 rp= RG.patchbase.first;
896                 while(rp) {
897                         next= rp->next;
898                         if(rp->f & RAD_SUBDIV) {
899                                 if(rp->emit[0]+rp->emit[1]+rp->emit[2]==0.0) {
900                                         contin= 1;
901                                         subdivideNode(rp->first, 0);
902                                         if(rp->first->down1) {
903                                                 subdivideNode(rp->first->down1, 0);
904                                                 subdivideNode(rp->first->down2, 0);
905                                         }
906                                 }
907                         }
908                         rp= next;
909                 }
910                 
911                 converttopatches();
912                 makeGlobalElemArray();
913
914                 if(contin==0 || MISC_test_break()) break;
915         }
916         make_node_display();
917 }
918
919 void inithemiwindows()
920 {
921         RadView *vw;
922
923         /* de hemiwindows */
924         vw= &(hemitop);
925         memset(vw, 0, sizeof(RadView));
926         vw->rectx= RG.hemires;
927         vw->recty= RG.hemires;
928         vw->rectz= MEM_mallocN(4*vw->rectx*vw->recty, "initwindows");
929         vw->rect= MEM_mallocN(4*vw->rectx*vw->recty, "initwindows");
930         vw->mynear= RG.maxsize/2000.0;
931         vw->myfar= 2.0*RG.maxsize;
932         vw->wx1= -vw->mynear;
933         vw->wx2= vw->mynear;
934         vw->wy1= -vw->mynear;
935         vw->wy2= vw->mynear;
936         
937         i_window(vw->wx1, vw->wx2, vw->wy1, vw->wy2, vw->mynear, vw->myfar, vw->winmat);
938
939         hemiside= hemitop;
940
941         vw= &(hemiside);
942         vw->recty/= 2;
943         vw->wy1= vw->wy2;
944         vw->wy2= 0.0;
945
946         i_window(vw->wx1, vw->wx2, vw->wy1, vw->wy2, vw->mynear, vw->myfar, vw->winmat);
947
948 }
949
950 void closehemiwindows()
951 {
952
953         if(hemiside.rect) MEM_freeN(hemiside.rect);
954         if(hemiside.rectz) MEM_freeN(hemiside.rectz);
955         hemiside.rectz= 0;
956         hemiside.rect= 0;
957         hemitop.rectz= 0;
958         hemitop.rect= 0;
959 }