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