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