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