resolved conflict state with HEAD r14096
[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                 Normalize(vec[0]);
86                 Normalize(vec[1]);
87                 Normalize(vec[2]);
88                 Normalize(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                 Normalize(cross[0]);
96                 Normalize(cross[1]);
97                 Normalize(cross[2]);
98                 Normalize(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                 Normalize(vec[0]);
130                 Normalize(vec[1]);
131                 Normalize(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                 Normalize(cross[0]);
138                 Normalize(cross[1]);
139                 Normalize(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= Normalize(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                                         INIT_MINMAX(min, max);
748                                         /* errmin and max are the filtered colors */
749                                         INIT_MINMAX(errmin, errmax);
750                                         minmaxradelemfilt(rp->first, min, max, errmin, errmax);
751                                         
752                                         /* if small difference between colors: no subdiv */
753                                         /* also test for the filtered ones: but with higher critical level */
754                                         
755                                         contin= 0;
756                                         a= abs( calculatecolor(min[0])-calculatecolor(max[0]));
757                                         b= abs( calculatecolor(errmin[0])-calculatecolor(errmax[0]));
758                                         if(a<15 || b<7) {
759                                                 c= abs( calculatecolor(min[1])-calculatecolor(max[1]));
760                                                 d= abs( calculatecolor(errmin[1])-calculatecolor(errmax[1]));
761                                                 if(c<15 || d<7) {
762                                                         e= abs( calculatecolor(min[2])-calculatecolor(max[2]));
763                                                         f= abs( calculatecolor(errmin[2])-calculatecolor(errmax[2]));
764                                                         if(e<15 || f<7) {
765                                                                 contin= 1;
766                                                                 clearsubflagelem(rp->first);
767                                                                 /* printf("%d %d %d %d %d %d\n", a, b, c, d, e, f); */
768                                                         }
769                                                 }
770                                         }
771                                         if(contin) {
772                                                 drawpatch_ext(rp, 0xFFFF);
773                                         }
774                                 }
775                         }
776                         rp->f &= ~RAD_SUBDIV;
777                         rp= rp->next;
778                 }
779                 
780                 contin= 0;
781                 
782                 el= RG.elem;
783                 for(b=RG.totelem; b>0; b--, el++) {
784                         rn= *el;
785                         if(rn->f & RAD_SUBDIV) {
786                                 rn->f-= RAD_SUBDIV;
787                                 subdivideNode(rn, 0);
788                                 if(rn->down1) {
789                                         subdivideNode(rn->down1, 0);
790                                         subdivideNode(rn->down2, 0);
791                                         contin= 1;
792                                 }
793                         }
794                 }
795                 makeGlobalElemArray();
796
797                 if(contin==0 || blender_test_break()) break;
798         }
799         
800         make_node_display();
801 }
802
803 void subdivideshootPatches(int it)
804 {
805         RPatch *rp, *shoot, *next;
806         float *fp, err, stoke, area;
807         int a, contin;
808         int maxlamp;
809         
810         if(RG.maxsublamp==0) maxlamp= RG.totlamp;
811         else maxlamp= RG.maxsublamp;
812         
813         while(it) {
814                 rad_printstatus();
815                 rad_init_energy();
816                 it--;
817                 
818                 for(a=0; a<maxlamp; a++) {
819                         shoot= findshootpatch();
820                         if(shoot==0) break;
821                         
822                         set_timecursor(a);
823                         drawpatch_ext(shoot, 0x88FF00);
824                         
825                         setnodeflags(shoot->first, RAD_SHOOT, 1);
826                         
827                         if( makeformfactors(shoot) ) {
828                         
829                                 fp= RG.formfactors;
830                                 rp= RG.patchbase.first;
831                                 while(rp) {
832                                         if(*fp!=0.0 && rp!=shoot) {
833                                         
834                                                 stoke= calcStokefactor(shoot, rp, rp->first, &area);
835                                                 if(stoke!= 0.0) {
836                                                         if(area>.1) {   /* does patch receive more than (about)10% of energy? */
837                                                                 rp->f= RAD_SUBDIV;
838                                                         }
839                                                         else {
840                                         
841                                                                 err= *fp/stoke;
842                                                         
843                                                                 /* area error */
844                                                                 area*=(0.5*RG.hemires*RG.hemires);
845                                                                 
846                                                                 if(area>45.0) {
847                                                                         if(err<0.95 || err>1.05) {
848                                                                                 if(err>0.05) {
849                                                                                                                 
850                                                                                         rp->f= RAD_SUBDIV;
851                                                                                         
852                                                                                 }
853                                                                         }
854                                                                 }
855                                                         }
856                                                 }
857                                         }
858                                         fp++;
859         
860                                         rp= rp->next;
861                                 }
862         
863                                 applyformfactors(shoot);
864                 
865                                 if( (a & 3)==1 ) {
866                                         make_node_display();
867                                         rad_forcedraw();
868                                 }
869                                 
870                                 setnodeflags(shoot->first, RAD_SHOOT, 0);
871                                 
872                                 if(blender_test_break()) break;
873                         }
874                         else a--;
875                         
876                 }
877                 
878                 contin= 0;
879
880                 rp= RG.patchbase.first;
881                 while(rp) {
882                         next= rp->next;
883                         if(rp->f & RAD_SUBDIV) {
884                                 if(rp->emit[0]+rp->emit[1]+rp->emit[2]==0.0) {
885                                         contin= 1;
886                                         subdivideNode(rp->first, 0);
887                                         if(rp->first->down1) {
888                                                 subdivideNode(rp->first->down1, 0);
889                                                 subdivideNode(rp->first->down2, 0);
890                                         }
891                                 }
892                         }
893                         rp= next;
894                 }
895                 
896                 converttopatches();
897                 makeGlobalElemArray();
898
899                 if(contin==0 || blender_test_break()) break;
900         }
901         make_node_display();
902 }
903
904 void inithemiwindows()
905 {
906         RadView *vw;
907
908         /* the hemiwindows */
909         vw= &(hemitop);
910         memset(vw, 0, sizeof(RadView));
911         vw->rectx= RG.hemires;
912         vw->recty= RG.hemires;
913         vw->rectz= MEM_mallocN(sizeof(int)*vw->rectx*vw->recty, "initwindows");
914         vw->rect= MEM_mallocN(sizeof(int)*vw->rectx*vw->recty, "initwindows");
915         vw->mynear= RG.maxsize/2000.0;
916         vw->myfar= 2.0*RG.maxsize;
917         vw->wx1= -vw->mynear;
918         vw->wx2= vw->mynear;
919         vw->wy1= -vw->mynear;
920         vw->wy2= vw->mynear;
921         
922         i_window(vw->wx1, vw->wx2, vw->wy1, vw->wy2, vw->mynear, vw->myfar, vw->winmat);
923
924         hemiside= hemitop;
925
926         vw= &(hemiside);
927         vw->recty/= 2;
928         vw->wy1= vw->wy2;
929         vw->wy2= 0.0;
930
931         i_window(vw->wx1, vw->wx2, vw->wy1, vw->wy2, vw->mynear, vw->myfar, vw->winmat);
932
933 }
934
935 void closehemiwindows()
936 {
937
938         if(hemiside.rect) MEM_freeN(hemiside.rect);
939         if(hemiside.rectz) MEM_freeN(hemiside.rectz);
940         hemiside.rectz= 0;
941         hemiside.rect= 0;
942         hemitop.rectz= 0;
943         hemitop.rect= 0;
944 }