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