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