Deep Shadow Buffer
[blender.git] / source / blender / render / intern / source / shadbuf.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * Contributor(s): 2004-2006, Blender Foundation
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  */
25
26 #include <math.h>
27 #include <string.h>
28
29
30 #include "MEM_guardedalloc.h"
31
32 #include "DNA_group_types.h"
33 #include "DNA_lamp_types.h"
34 #include "DNA_material_types.h"
35
36 #include "BKE_global.h"
37 #include "BKE_scene.h"
38 #include "BKE_utildefines.h"
39
40 #include "BLI_arithb.h"
41 #include "BLI_blenlib.h"
42 #include "BLI_jitter.h"
43 #include "BLI_memarena.h"
44 #include "BLI_rand.h"
45
46 #include "PIL_time.h"
47
48 #include "renderpipeline.h"
49 #include "render_types.h"
50 #include "renderdatabase.h"
51 #include "rendercore.h"
52 #include "shadbuf.h"
53 #include "shading.h"
54 #include "zbuf.h"
55
56 /* XXX, could be better implemented... this is for endian issues
57 */
58 #if defined(__sgi) || defined(__sparc) || defined(__sparc__) || defined (__PPC__) || defined (__ppc__) || defined (__hppa__) || defined (__BIG_ENDIAN__)
59 #define RCOMP   3
60 #define GCOMP   2
61 #define BCOMP   1
62 #define ACOMP   0
63 #else
64 #define RCOMP   0
65 #define GCOMP   1
66 #define BCOMP   2
67 #define ACOMP   3
68 #endif
69
70 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
71 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
72 /* only to be used here in this file, it's for speed */
73 extern struct Render R;
74 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
75
76 /* ------------------------------------------------------------------------- */
77
78 /* initshadowbuf() in convertBlenderScene.c */
79
80 /* ------------------------------------------------------------------------- */
81
82 static void copy_to_ztile(int *rectz, int size, int x1, int y1, int tile, char *r1)
83 {
84         int len4, *rz;  
85         int x2, y2;
86         
87         x2= x1+tile;
88         y2= y1+tile;
89         if(x2>=size) x2= size-1;
90         if(y2>=size) y2= size-1;
91
92         if(x1>=x2 || y1>=y2) return;
93
94         len4= 4*(x2- x1);
95         rz= rectz + size*y1 + x1;
96         for(; y1<y2; y1++) {
97                 memcpy(r1, rz, len4);
98                 rz+= size;
99                 r1+= len4;
100         }
101 }
102
103 #if 0
104 static int sizeoflampbuf(ShadBuf *shb)
105 {
106         int num,count=0;
107         char *cp;
108         
109         cp= shb->cbuf;
110         num= (shb->size*shb->size)/256;
111
112         while(num--) count+= *(cp++);
113         
114         return 256*count;
115 }
116 #endif
117
118 /* not threadsafe... */
119 static float *give_jitter_tab(int samp)
120 {
121         /* these are all possible jitter tables, takes up some
122          * 12k, not really bad!
123          * For soft shadows, it saves memory and render time
124          */
125         static int tab[17]={1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, 196, 225, 256};
126         static float jit[1496][2];
127         static char ctab[17]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
128         int a, offset=0;
129         
130         if(samp<2) samp= 2;
131         else if(samp>16) samp= 16;
132
133         for(a=0; a<samp-1; a++) offset+= tab[a];
134
135         if(ctab[samp]==0) {
136                 ctab[samp]= 1;
137                 BLI_initjit(jit[offset], samp*samp);
138         }
139                 
140         return jit[offset];
141         
142 }
143
144 static void make_jitter_weight_tab(Render *re, ShadBuf *shb, short filtertype) 
145 {
146         float *jit, totw= 0.0f;
147         int samp= get_render_shadow_samples(&re->r, shb->samp);
148         int a, tot=samp*samp;
149         
150         shb->weight= MEM_mallocN(sizeof(float)*tot, "weight tab lamp");
151         
152         for(jit= shb->jit, a=0; a<tot; a++, jit+=2) {
153                 if(filtertype==LA_SHADBUF_TENT) 
154                         shb->weight[a]= 0.71f - sqrt(jit[0]*jit[0] + jit[1]*jit[1]);
155                 else if(filtertype==LA_SHADBUF_GAUSS) 
156                         shb->weight[a]= RE_filter_value(R_FILTER_GAUSS, 1.8f*sqrt(jit[0]*jit[0] + jit[1]*jit[1]));
157                 else
158                         shb->weight[a]= 1.0f;
159                 
160                 totw+= shb->weight[a];
161         }
162         
163         totw= 1.0f/totw;
164         for(a=0; a<tot; a++) {
165                 shb->weight[a]*= totw;
166         }
167 }
168
169 static int verg_deepsample(const void *poin1, const void *poin2)
170 {
171         const DeepSample *ds1= (const DeepSample*)poin1;
172         const DeepSample *ds2= (const DeepSample*)poin2;
173
174         if(ds1->z < ds2->z) return -1;
175         else if(ds1->z == ds2->z) return 0;
176         else return 1;
177 }
178
179 static int compress_deepsamples(DeepSample *dsample, int tot, float epsilon)
180 {
181         /* uses doubles to avoid overflows and other numerical issues,
182            could be improved */
183         DeepSample *ds, *newds;
184         float v;
185         double slope, slopemin, slopemax, min, max, div, newmin, newmax;
186         int a, first, z, newtot= 0;
187
188         /*if(print) {
189                 for(a=0, ds=dsample; a<tot; a++, ds++)
190                         printf("%lf,%f ", ds->z/(double)0x7FFFFFFF, ds->v);
191                 printf("\n");
192         }*/
193
194         /* read from and write into same array */
195         ds= dsample;
196         newds= dsample;
197         a= 0;
198
199         /* as long as we are not at the end of the array */
200         for(a++, ds++; a<tot; a++, ds++) {
201                 slopemin= 0.0f;
202                 slopemax= 0.0f;
203                 first= 1;
204
205                 for(; a<tot; a++, ds++) {
206                         //dz= ds->z - newds->z;
207                         if(ds->z == newds->z) {
208                                 /* still in same z position, simply check
209                                    visibility difference against epsilon */
210                                 if(!(fabs(newds->v - ds->v) <= epsilon)) {
211                                         break;
212                                 }
213                         }
214                         else {
215                                 /* compute slopes */
216                                 div= (double)0x7FFFFFFF/((double)ds->z - (double)newds->z);
217                                 min= ((ds->v - epsilon) - newds->v)*div;
218                                 max= ((ds->v + epsilon) - newds->v)*div;
219
220                                 /* adapt existing slopes */
221                                 if(first) {
222                                         newmin= min;
223                                         newmax= max;
224                                         first= 0;
225                                 }
226                                 else {
227                                         newmin= MAX2(slopemin, min);
228                                         newmax= MIN2(slopemax, max);
229
230                                         /* verify if there is still space between the slopes */
231                                         if(newmin > newmax) {
232                                                 ds--;
233                                                 a--;
234                                                 break;
235                                         }
236                                 }
237
238                                 slopemin= newmin;
239                                 slopemax= newmax;
240                         }
241                 }
242
243                 if(a == tot) {
244                         ds--;
245                         a--;
246                 }
247
248                 /* always previous z */
249                 z= ds->z;
250
251                 if(first || a==tot-1) {
252                         /* if slopes were not initialized, use last visibility */
253                         v= ds->v;
254                 }
255                 else {
256                         /* compute visibility at center between slopes at z */
257                         slope= (slopemin+slopemax)*0.5;
258                         v= newds->v + slope*((z - newds->z)/(double)0x7FFFFFFF);
259                 }
260
261                 newds++;
262                 newtot++;
263
264                 newds->z= z;
265                 newds->v= v;
266         }
267
268         if(newtot == 0 || (newds->v != (newds-1)->v))
269                 newtot++;
270
271         /*if(print) {
272                 for(a=0, ds=dsample; a<newtot; a++, ds++)
273                         printf("%lf,%f ", ds->z/(double)0x7FFFFFFF, ds->v);
274                 printf("\n");
275         }*/
276
277         return newtot;
278 }
279
280 static float deep_alpha(Render *re, int obinr, int facenr, int strand)
281 {
282         ObjectInstanceRen *obi= &re->objectinstance[obinr];
283         Material *ma;
284
285         if(strand) {
286                 StrandRen *strand= RE_findOrAddStrand(obi->obr, facenr-1);
287                 ma= strand->buffer->ma;
288         }
289         else {
290                 VlakRen *vlr= RE_findOrAddVlak(obi->obr, (facenr-1) & RE_QUAD_MASK);
291                 ma= vlr->mat;
292         }
293
294         return ma->shad_alpha;
295 }
296
297 static void compress_deepshadowbuf(Render *re, ShadBuf *shb, APixstr *apixbuf, APixstrand *apixbufstrand)
298 {
299         ShadSampleBuf *shsample;
300         DeepSample *ds[RE_MAX_OSA], *sampleds[RE_MAX_OSA], *dsb, *newbuf;
301         APixstr *ap, *apn;
302         APixstrand *aps, *apns;
303         float visibility, totbuf= shb->totbuf;
304         int a, b, c, tot, minz, found, size= shb->size, prevtot, newtot;
305         int sampletot[RE_MAX_OSA], totsample = 0, totsamplec = 0;
306         
307         shsample= MEM_callocN( sizeof(ShadSampleBuf), "shad sample buf");
308         BLI_addtail(&shb->buffers, shsample);
309
310         shsample->totbuf= MEM_callocN(sizeof(int)*size*size, "deeptotbuf");
311         shsample->deepbuf= MEM_callocN(sizeof(DeepSample*)*size*size, "deepbuf");
312
313         ap= apixbuf;
314         aps= apixbufstrand;
315         for(a=0; a<size*size; a++, ap++, aps++) {
316                 /* count number of samples */
317                 for(c=0; c<totbuf; c++)
318                         sampletot[c]= 0;
319
320                 tot= 0;
321                 for(apn=ap; apn; apn=apn->next)
322                         for(b=0; b<4; b++)
323                                 if(apn->p[b])
324                                         for(c=0; c<totbuf; c++)
325                                                 if(apn->mask[b] & (1<<c))
326                                                         sampletot[c]++;
327
328                 if(apixbufstrand) {
329                         for(apns=aps; apns; apns=apns->next)
330                                 for(b=0; b<4; b++)
331                                         if(apns->p[b])
332                                                 for(c=0; c<totbuf; c++)
333                                                         if(apns->mask[b] & (1<<c))
334                                                                 sampletot[c]++;
335                 }
336
337                 for(c=0; c<totbuf; c++)
338                         tot += sampletot[c];
339
340                 if(tot == 0) {
341                         shsample->deepbuf[a]= NULL;
342                         shsample->totbuf[a]= 0;
343                         continue;
344                 }
345
346                 /* fill samples */
347                 ds[0]= sampleds[0]= MEM_callocN(sizeof(DeepSample)*tot*2, "deepsample");
348                 for(c=1; c<totbuf; c++)
349                         ds[c]= sampleds[c]= sampleds[c-1] + sampletot[c-1]*2;
350
351                 for(apn=ap; apn; apn=apn->next) {
352                         for(b=0; b<4; b++) {
353                                 if(apn->p[b]) {
354                                         for(c=0; c<totbuf; c++) {
355                                                 if(apn->mask[b] & (1<<c)) {
356                                                         /* two entries to create step profile */
357                                                         ds[c]->z= apn->z[b];
358                                                         ds[c]->v= 1.0f; /* not used */
359                                                         ds[c]++;
360                                                         ds[c]->z= apn->z[b];
361                                                         ds[c]->v= deep_alpha(re, apn->obi[b], apn->p[b], 0);
362                                                         ds[c]++;
363                                                 }
364                                         }
365                                 }
366                         }
367                 }
368
369                 if(apixbufstrand) {
370                         for(apns=aps; apns; apns=apns->next) {
371                                 for(b=0; b<4; b++) {
372                                         if(apns->p[b]) {
373                                                 for(c=0; c<totbuf; c++) {
374                                                         if(apns->mask[b] & (1<<c)) {
375                                                                 /* two entries to create step profile */
376                                                                 ds[c]->z= apns->z[b];
377                                                                 ds[c]->v= 1.0f; /* not used */
378                                                                 ds[c]++;
379                                                                 ds[c]->z= apns->z[b];
380                                                                 ds[c]->v= deep_alpha(re, apns->obi[b], apns->p[b], 1);
381                                                                 ds[c]++;
382                                                         }
383                                                 }
384                                         }
385                                 }
386                         }
387                 }
388
389                 for(c=0; c<totbuf; c++) {
390                         /* sort by increasing z */
391                         qsort(sampleds[c], sampletot[c], sizeof(DeepSample)*2, verg_deepsample);
392
393                         /* sum visibility, replacing alpha values */
394                         visibility= 1.0f;
395                         ds[c]= sampleds[c];
396
397                         for(b=0; b<sampletot[c]; b++) {
398                                 /* two entries creating step profile */
399                                 ds[c]->v= visibility;
400                                 ds[c]++;
401
402                                 visibility *= 1.0f-ds[c]->v;
403                                 ds[c]->v= visibility;
404                                 ds[c]++;
405                         }
406
407                         /* halfway trick, probably won't work well for volumes? */
408                         ds[c]= sampleds[c];
409                         for(b=0; b<sampletot[c]; b++) {
410                                 if(b+1 < sampletot[c]) {
411                                         ds[c]->z= (ds[c]->z>>1) + ((ds[c]+2)->z>>1);
412                                         ds[c]++;
413                                         ds[c]->z= (ds[c]->z>>1) + ((ds[c]+2)->z>>1);
414                                         ds[c]++;
415                                 }
416                                 else {
417                                         ds[c]->z= (ds[c]->z>>1) + (0x7FFFFFFF>>1);
418                                         ds[c]++;
419                                         ds[c]->z= (ds[c]->z>>1) + (0x7FFFFFFF>>1);
420                                         ds[c]++;
421                                 }
422                         }
423
424                         /* init for merge loop */
425                         ds[c]= sampleds[c];
426                         sampletot[c] *= 2;
427                 }
428
429                 shsample->deepbuf[a]= MEM_callocN(sizeof(DeepSample)*tot*2, "deepsample");
430                 shsample->totbuf[a]= 0;
431
432                 /* merge buffers */
433                 dsb= shsample->deepbuf[a];
434                 while(1) {
435                         minz= 0;
436                         found= 0;
437
438                         for(c=0; c<totbuf; c++) {
439                                 if(sampletot[c] && (!found || ds[c]->z < minz)) {
440                                         minz= ds[c]->z;
441                                         found= 1;
442                                 }
443                         }
444
445                         if(!found)
446                                 break;
447
448                         dsb->z= minz;
449                         dsb->v= 0.0f;
450
451                         visibility= 0.0f;
452                         for(c=0; c<totbuf; c++) {
453                                 if(sampletot[c] && ds[c]->z == minz) {
454                                         ds[c]++;
455                                         sampletot[c]--;
456                                 }
457
458                                 if(sampleds[c] == ds[c])
459                                         visibility += 1.0f/totbuf;
460                                 else
461                                         visibility += (ds[c]-1)->v/totbuf;
462                         }
463
464                         dsb->v= visibility;
465                         dsb++;
466                         shsample->totbuf[a]++;
467                 }
468
469                 prevtot= shsample->totbuf[a];
470                 totsample += prevtot;
471
472                 newtot= compress_deepsamples(shsample->deepbuf[a], prevtot, shb->compressthresh);
473                 shsample->totbuf[a]= newtot;
474                 totsamplec += newtot;
475
476                 if(newtot < prevtot) {
477                         newbuf= MEM_mallocN(sizeof(DeepSample)*newtot, "cdeepsample");
478                         memcpy(newbuf, shsample->deepbuf[a], sizeof(DeepSample)*newtot);
479                         MEM_freeN(shsample->deepbuf[a]);
480                         shsample->deepbuf[a]= newbuf;
481                 }
482
483                 MEM_freeN(sampleds[0]);
484         }
485
486         //printf("%d -> %d, ratio %f\n", totsample, totsamplec, (float)totsamplec/(float)totsample);
487 }
488
489 /* create Z tiles (for compression): this system is 24 bits!!! */
490 static void compress_shadowbuf(ShadBuf *shb, int *rectz, int square)
491 {
492         ShadSampleBuf *shsample;
493         float dist;
494         uintptr_t *ztile;
495         int *rz, *rz1, verg, verg1, size= shb->size;
496         int a, x, y, minx, miny, byt1, byt2;
497         char *rc, *rcline, *ctile, *zt;
498         
499         shsample= MEM_callocN( sizeof(ShadSampleBuf), "shad sample buf");
500         BLI_addtail(&shb->buffers, shsample);
501         
502         shsample->zbuf= MEM_mallocN( sizeof(uintptr_t)*(size*size)/256, "initshadbuf2");
503         shsample->cbuf= MEM_callocN( (size*size)/256, "initshadbuf3");
504         
505         ztile= (uintptr_t *)shsample->zbuf;
506         ctile= shsample->cbuf;
507         
508         /* help buffer */
509         rcline= MEM_mallocN(256*4+sizeof(int), "makeshadbuf2");
510         
511         for(y=0; y<size; y+=16) {
512                 if(y< size/2) miny= y+15-size/2;
513                 else miny= y-size/2;    
514                 
515                 for(x=0; x<size; x+=16) {
516                         
517                         /* is tile within spotbundle? */
518                         a= size/2;
519                         if(x< a) minx= x+15-a;
520                         else minx= x-a; 
521                         
522                         dist= sqrt( (float)(minx*minx+miny*miny) );
523                         
524                         if(square==0 && dist>(float)(a+12)) {   /* 12, tested with a onlyshadow lamp */
525                                 a= 256; verg= 0; /* 0x80000000; */ /* 0x7FFFFFFF; */
526                                 rz1= (&verg)+1;
527                         } 
528                         else {
529                                 copy_to_ztile(rectz, size, x, y, 16, rcline);
530                                 rz1= (int *)rcline;
531                                 
532                                 verg= (*rz1 & 0xFFFFFF00);
533                                 
534                                 for(a=0;a<256;a++,rz1++) {
535                                         if( (*rz1 & 0xFFFFFF00) !=verg) break;
536                                 }
537                         }
538                         if(a==256) { /* complete empty tile */
539                                 *ctile= 0;
540                                 *ztile= *(rz1-1);
541                         }
542                         else {
543                                 
544                                 /* ACOMP etc. are defined to work L/B endian */
545                                 
546                                 rc= rcline;
547                                 rz1= (int *)rcline;
548                                 verg=  rc[ACOMP];
549                                 verg1= rc[BCOMP];
550                                 rc+= 4;
551                                 byt1= 1; byt2= 1;
552                                 for(a=1;a<256;a++,rc+=4) {
553                                         byt1 &= (verg==rc[ACOMP]);
554                                         byt2 &= (verg1==rc[BCOMP]);
555                                         
556                                         if(byt1==0) break;
557                                 }
558                                 if(byt1 && byt2) {      /* only store byte */
559                                         *ctile= 1;
560                                         *ztile= (uintptr_t)MEM_mallocN(256+4, "tile1");
561                                         rz= (int *)*ztile;
562                                         *rz= *rz1;
563                                         
564                                         zt= (char *)(rz+1);
565                                         rc= rcline;
566                                         for(a=0; a<256; a++, zt++, rc+=4) *zt= rc[GCOMP];       
567                                 }
568                                 else if(byt1) {         /* only store short */
569                                         *ctile= 2;
570                                         *ztile= (uintptr_t)MEM_mallocN(2*256+4,"Tile2");
571                                         rz= (int *)*ztile;
572                                         *rz= *rz1;
573                                         
574                                         zt= (char *)(rz+1);
575                                         rc= rcline;
576                                         for(a=0; a<256; a++, zt+=2, rc+=4) {
577                                                 zt[0]= rc[BCOMP];
578                                                 zt[1]= rc[GCOMP];
579                                         }
580                                 }
581                                 else {                  /* store triple */
582                                         *ctile= 3;
583                                         *ztile= (uintptr_t)MEM_mallocN(3*256,"Tile3");
584
585                                         zt= (char *)*ztile;
586                                         rc= rcline;
587                                         for(a=0; a<256; a++, zt+=3, rc+=4) {
588                                                 zt[0]= rc[ACOMP];
589                                                 zt[1]= rc[BCOMP];
590                                                 zt[2]= rc[GCOMP];
591                                         }
592                                 }
593                         }
594                         ztile++;
595                         ctile++;
596                 }
597         }
598
599         MEM_freeN(rcline);
600 }
601
602 /* sets start/end clipping. lar->shb should be initialized */
603 static void shadowbuf_autoclip(Render *re, LampRen *lar)
604 {
605         ObjectInstanceRen *obi;
606         ObjectRen *obr;
607         VlakRen *vlr= NULL;
608         VertRen *ver= NULL;
609         Material *ma= NULL;
610         float minz, maxz, vec[3], viewmat[4][4], obviewmat[4][4];
611         unsigned int lay = -1;
612         int i, a, maxtotvert, ok= 1;
613         char *clipflag;
614         
615         minz= 1.0e30f; maxz= -1.0e30f;
616         Mat4CpyMat4(viewmat, lar->shb->viewmat);
617         
618         if(lar->mode & (LA_LAYER|LA_LAYER_SHADOW)) lay= lar->lay;
619
620         maxtotvert= 0;
621         for(obr=re->objecttable.first; obr; obr=obr->next)
622                 maxtotvert= MAX2(obr->totvert, maxtotvert);
623
624         clipflag= MEM_callocN(sizeof(char)*maxtotvert, "autoclipflag");
625
626         /* set clip in vertices when face visible */
627         for(i=0, obi=re->instancetable.first; obi; i++, obi=obi->next) {
628                 obr= obi->obr;
629
630                 if(obi->flag & R_TRANSFORMED)
631                         Mat4MulMat4(obviewmat, obi->mat, viewmat);
632                 else
633                         Mat4CpyMat4(obviewmat, viewmat);
634
635                 memset(clipflag, 0, sizeof(char)*obr->totvert);
636
637                 /* clear clip, is being set if face is visible (clip is calculated for real later) */
638                 for(a=0; a<obr->totvlak; a++) {
639                         if((a & 255)==0) vlr= obr->vlaknodes[a>>8].vlak;
640                         else vlr++;
641                         
642                         /* note; these conditions are copied from zbuffer_shadow() */
643                         if(vlr->mat!= ma) {
644                                 ma= vlr->mat;
645                                 ok= 1;
646                                 if((ma->mode & MA_SHADBUF)==0) ok= 0;
647                         }
648                         
649                         if(ok && (obi->lay & lay)) {
650                                 clipflag[vlr->v1->index]= 1;
651                                 clipflag[vlr->v2->index]= 1;
652                                 clipflag[vlr->v3->index]= 1;
653                                 if(vlr->v4) clipflag[vlr->v4->index]= 1;
654                         }                               
655                 }               
656                 
657                 /* calculate min and max */
658                 for(a=0; a< obr->totvert;a++) {
659                         if((a & 255)==0) ver= RE_findOrAddVert(obr, a);
660                         else ver++;
661                         
662                         if(clipflag[a]) {
663                                 VECCOPY(vec, ver->co);
664                                 Mat4MulVecfl(obviewmat, vec);
665                                 /* Z on visible side of lamp space */
666                                 if(vec[2] < 0.0f) {
667                                         float inpr, z= -vec[2];
668                                         
669                                         /* since vec is rotated in lampspace, this is how to get the cosine of angle */
670                                         /* precision is set 20% larger */
671                                         vec[2]*= 1.2f;
672                                         Normalize(vec);
673                                         inpr= - vec[2];
674
675                                         if(inpr>=lar->spotsi) {
676                                                 if(z<minz) minz= z;
677                                                 if(z>maxz) maxz= z;
678                                         }
679                                 }
680                         }
681                 }
682         }
683
684         MEM_freeN(clipflag);
685         
686         /* set clipping min and max */
687         if(minz < maxz) {
688                 float delta= (maxz - minz);     /* threshold to prevent precision issues */
689                 
690                 //printf("minz %f maxz %f delta %f\n", minz, maxz, delta);
691                 if(lar->bufflag & LA_SHADBUF_AUTO_START)
692                         lar->shb->d= minz - delta*0.02f;        /* 0.02 is arbitrary... needs more thinking! */
693                 if(lar->bufflag & LA_SHADBUF_AUTO_END)
694                         lar->shb->clipend= maxz + delta*0.1f;
695                 
696                 /* bias was calculated as percentage, we scale it to prevent animation issues */
697                 delta= (lar->clipend-lar->clipsta)/(lar->shb->clipend-lar->shb->d);
698                 //printf("bias delta %f\n", delta);
699                 lar->shb->bias= (int) (delta*(float)lar->shb->bias);
700         }
701 }
702
703 static void makeflatshadowbuf(Render *re, LampRen *lar, float *jitbuf)
704 {
705         ShadBuf *shb= lar->shb;
706         int *rectz, samples;
707
708         /* zbuffering */
709         rectz= MEM_mapallocN(sizeof(int)*shb->size*shb->size, "makeshadbuf");
710         
711         for(samples=0; samples<shb->totbuf; samples++) {
712                 zbuffer_shadow(re, shb->persmat, lar, rectz, shb->size, jitbuf[2*samples], jitbuf[2*samples+1]);
713                 /* create Z tiles (for compression): this system is 24 bits!!! */
714                 compress_shadowbuf(shb, rectz, lar->mode & LA_SQUARE);
715
716                 if(re->test_break(re->tbh))
717                         break;
718         }
719         
720         MEM_freeN(rectz);
721 }
722
723 static void makedeepshadowbuf(Render *re, LampRen *lar, float *jitbuf)
724 {
725         ShadBuf *shb= lar->shb;
726         APixstr *apixbuf;
727         APixstrand *apixbufstrand= NULL;
728         ListBase apsmbase= {NULL, NULL};
729
730         /* zbuffering */
731         apixbuf= MEM_callocN(sizeof(APixstr)*shb->size*shb->size, "APixbuf");
732         if(re->totstrand)
733                 apixbufstrand= MEM_callocN(sizeof(APixstrand)*shb->size*shb->size, "APixbufstrand");
734
735         zbuffer_abuf_shadow(re, lar, shb->persmat, apixbuf, apixbufstrand, &apsmbase, shb->size,
736                 shb->totbuf, (float(*)[2])jitbuf);
737
738         /* create Z tiles (for compression): this system is 24 bits!!! */
739         compress_deepshadowbuf(re, shb, apixbuf, apixbufstrand);
740         
741         MEM_freeN(apixbuf);
742         if(apixbufstrand)
743                 MEM_freeN(apixbufstrand);
744         freepsA(&apsmbase);
745 }
746
747 void makeshadowbuf(Render *re, LampRen *lar)
748 {
749         ShadBuf *shb= lar->shb;
750         float wsize, *jitbuf, twozero[2]= {0.0f, 0.0f}, angle, temp;
751         
752         if(lar->bufflag & (LA_SHADBUF_AUTO_START|LA_SHADBUF_AUTO_END))
753                 shadowbuf_autoclip(re, lar);
754         
755         /* just to enforce identical behaviour of all irregular buffers */
756         if(lar->buftype==LA_SHADBUF_IRREGULAR)
757                 shb->size= 1024;
758         
759         /* matrices and window: in winmat the transformation is being put,
760                 transforming from observer view to lamp view, including lamp window matrix */
761         
762         angle= saacos(lar->spotsi);
763         temp= 0.5f*shb->size*cos(angle)/sin(angle);
764         shb->pixsize= (shb->d)/temp;
765         wsize= shb->pixsize*(shb->size/2.0);
766         
767         i_window(-wsize, wsize, -wsize, wsize, shb->d, shb->clipend, shb->winmat);
768         Mat4MulMat4(shb->persmat, shb->viewmat, shb->winmat);
769
770         if(ELEM3(lar->buftype, LA_SHADBUF_REGULAR, LA_SHADBUF_HALFWAY, LA_SHADBUF_DEEP)) {
771                 shb->totbuf= lar->buffers;
772
773                 /* jitter, weights - not threadsafe! */
774                 BLI_lock_thread(LOCK_CUSTOM1);
775                 shb->jit= give_jitter_tab(get_render_shadow_samples(&re->r, shb->samp));
776                 make_jitter_weight_tab(re, shb, lar->filtertype);
777                 BLI_unlock_thread(LOCK_CUSTOM1);
778                 
779                 if(shb->totbuf==4) jitbuf= give_jitter_tab(2);
780                 else if(shb->totbuf==9) jitbuf= give_jitter_tab(3);
781                 else jitbuf= twozero;
782                 
783                 /* zbuffering */
784                 if(lar->buftype == LA_SHADBUF_DEEP) {
785                         makedeepshadowbuf(re, lar, jitbuf);
786                         shb->totbuf= 1;
787                 }
788                 else
789                         makeflatshadowbuf(re, lar, jitbuf);
790
791                 /* printf("lampbuf %d\n", sizeoflampbuf(shb)); */
792         }
793 }
794
795 static void *do_shadow_thread(void *re_v)
796 {
797         Render *re= (Render*)re_v;
798         LampRen *lar;
799
800         do {
801                 BLI_lock_thread(LOCK_CUSTOM1);
802                 for(lar=re->lampren.first; lar; lar=lar->next) {
803                         if(lar->shb && !lar->thread_assigned) {
804                                 lar->thread_assigned= 1;
805                                 break;
806                         }
807                 }
808                 BLI_unlock_thread(LOCK_CUSTOM1);
809
810                 /* if type is irregular, this only sets the perspective matrix and autoclips */
811                 if(lar) {
812                         makeshadowbuf(re, lar);
813                         BLI_lock_thread(LOCK_CUSTOM1);
814                         lar->thread_ready= 1;
815                         BLI_unlock_thread(LOCK_CUSTOM1);
816                 }
817         } while(lar && !re->test_break(re->tbh));
818
819         return NULL;
820 }
821
822 static volatile int g_break= 0;
823 static int thread_break(void *unused)
824 {
825         return g_break;
826 }
827
828 void threaded_makeshadowbufs(Render *re)
829 {
830         ListBase threads;
831         LampRen *lar;
832         int a, totthread= 0;
833         int (*test_break)(void *);
834
835         /* count number of threads to use */
836         if(G.rendering) {
837                 for(lar=re->lampren.first; lar; lar= lar->next)
838                         if(lar->shb)
839                                 totthread++;
840                 
841                 totthread= MIN2(totthread, re->r.threads);
842         }
843         else
844                 totthread= 1; /* preview render */
845
846         if(totthread <= 1) {
847                 for(lar=re->lampren.first; lar; lar= lar->next) {
848                         if(re->test_break(re->tbh)) break;
849                         if(lar->shb) {
850                                 /* if type is irregular, this only sets the perspective matrix and autoclips */
851                                 makeshadowbuf(re, lar);
852                         }
853                 }
854         }
855         else {
856                 /* swap test break function */
857                 test_break= re->test_break;
858                 re->test_break= thread_break;
859
860                 for(lar=re->lampren.first; lar; lar= lar->next) {
861                         lar->thread_assigned= 0;
862                         lar->thread_ready= 0;
863                 }
864
865                 BLI_init_threads(&threads, do_shadow_thread, totthread);
866                 
867                 for(a=0; a<totthread; a++)
868                         BLI_insert_thread(&threads, re);
869
870                 /* keep rendering as long as there are shadow buffers not ready */
871                 do {
872                         if((g_break=test_break(re->tbh)))
873                                 break;
874
875                         PIL_sleep_ms(50);
876
877                         BLI_lock_thread(LOCK_CUSTOM1);
878                         for(lar=re->lampren.first; lar; lar= lar->next)
879                                 if(lar->shb && !lar->thread_ready)
880                                         break;
881                         BLI_unlock_thread(LOCK_CUSTOM1);
882                 } while(lar);
883         
884                 BLI_end_threads(&threads);
885
886                 /* unset threadsafety */
887                 re->test_break= test_break;
888                 g_break= 0;
889         }
890 }
891
892 void freeshadowbuf(LampRen *lar)
893 {
894         if(lar->shb) {
895                 ShadBuf *shb= lar->shb;
896                 ShadSampleBuf *shsample;
897                 int b, v;
898                 
899                 for(shsample= shb->buffers.first; shsample; shsample= shsample->next) {
900                         if(shsample->deepbuf) {
901                                 v= shb->size*shb->size;
902                                 for(b=0; b<v; b++)
903                                         if(shsample->deepbuf[b])
904                                                 MEM_freeN(shsample->deepbuf[b]);
905                                         
906                                 MEM_freeN(shsample->deepbuf);
907                                 MEM_freeN(shsample->totbuf);
908                         }
909                         else {
910                                 intptr_t *ztile= shsample->zbuf;
911                                 char *ctile= shsample->cbuf;
912                                 
913                                 v= (shb->size*shb->size)/256;
914                                 for(b=0; b<v; b++, ztile++, ctile++)
915                                         if(*ctile) MEM_freeN((void *) *ztile);
916                                 
917                                 MEM_freeN(shsample->zbuf);
918                                 MEM_freeN(shsample->cbuf);
919                         }
920                 }
921                 BLI_freelistN(&shb->buffers);
922                 
923                 if(shb->weight) MEM_freeN(shb->weight);
924                 MEM_freeN(lar->shb);
925                 
926                 lar->shb= NULL;
927         }
928 }
929
930
931 static int firstreadshadbuf(ShadBuf *shb, ShadSampleBuf *shsample, int **rz, int xs, int ys, int nr)
932 {
933         /* return a 1 if fully compressed shadbuf-tile && z==const */
934         int ofs;
935         char *ct;
936
937         if(shsample->deepbuf)
938                 return 0;
939
940         /* always test borders of shadowbuffer */
941         if(xs<0) xs= 0; else if(xs>=shb->size) xs= shb->size-1;
942         if(ys<0) ys= 0; else if(ys>=shb->size) ys= shb->size-1;
943    
944         /* calc z */
945         ofs= (ys>>4)*(shb->size>>4) + (xs>>4);
946         ct= shsample->cbuf+ofs;
947         if(*ct==0) {
948             if(nr==0) {
949                         *rz= *( (int **)(shsample->zbuf+ofs) );
950                         return 1;
951             }
952                 else if(*rz!= *( (int **)(shsample->zbuf+ofs) )) return 0;
953                 
954             return 1;
955         }
956         
957         return 0;
958 }
959
960 static float readdeepvisibility(DeepSample *dsample, int tot, int z, int bias, float *biast)
961 {
962         DeepSample *ds, *prevds;
963         float t;
964         int a;
965
966         /* tricky stuff here; we use ints which can overflow easily with bias values */
967
968         ds= dsample;
969         for(a=0; a<tot && (z-bias > ds->z); a++, ds++)
970                 ;
971
972         if(a == tot) {
973                 if(biast)
974                         *biast= 0.0f;
975                 return (ds-1)->v; /* completely behind all samples */
976         }
977         
978         /* check if this read needs bias blending */
979         if(biast) {
980                 if(z > ds->z)
981                         *biast= (float)(z - ds->z)/(float)bias;
982                 else
983                         *biast= 0.0f;
984         }
985
986         if(a == 0)
987                 return 1.0f; /* completely in front of all samples */
988
989         prevds= ds-1;
990         t= (float)(z-bias - prevds->z)/(float)(ds->z - prevds->z);
991         return t*ds->v + (1.0f-t)*prevds->v;
992 }
993
994 static float readdeepshadowbuf(ShadBuf *shb, ShadSampleBuf *shsample, int bias, int xs, int ys, int zs)
995 {
996         float v, biasv, biast;
997         int ofs, tot;
998
999         if(zs < - 0x7FFFFE00 + bias)
1000                 return 1.0;     /* extreme close to clipstart */
1001
1002         /* calc z */
1003         ofs= ys*shb->size + xs;
1004         tot= shsample->totbuf[ofs];
1005         if(tot == 0)
1006                 return 1.0f;
1007
1008         v= readdeepvisibility(shsample->deepbuf[ofs], tot, zs, bias, &biast);
1009
1010         if(biast != 0.0f) {
1011                 /* in soft bias area */
1012                 biasv= readdeepvisibility(shsample->deepbuf[ofs], tot, zs, 0, 0);
1013
1014                 biast= biast*biast;
1015                 return (1.0f-biast)*v + biast*biasv;
1016         }
1017
1018         return v;
1019 }
1020
1021 /* return 1.0 : fully in light */
1022 static float readshadowbuf(ShadBuf *shb, ShadSampleBuf *shsample, int bias, int xs, int ys, int zs)     
1023 {
1024         float temp;
1025         int *rz, ofs;
1026         int zsamp=0;
1027         char *ct, *cz;
1028
1029         /* simpleclip */
1030         /* if(xs<0 || ys<0) return 1.0; */
1031         /* if(xs>=shb->size || ys>=shb->size) return 1.0; */
1032         
1033         /* always test borders of shadowbuffer */
1034         if(xs<0) xs= 0; else if(xs>=shb->size) xs= shb->size-1;
1035         if(ys<0) ys= 0; else if(ys>=shb->size) ys= shb->size-1;
1036
1037         if(shsample->deepbuf)
1038                 return readdeepshadowbuf(shb, shsample, bias, xs, ys, zs);
1039
1040         /* calc z */
1041         ofs= (ys>>4)*(shb->size>>4) + (xs>>4);
1042         ct= shsample->cbuf+ofs;
1043         rz= *( (int **)(shsample->zbuf+ofs) );
1044
1045         if(*ct==3) {
1046                 ct= ((char *)rz)+3*16*(ys & 15)+3*(xs & 15);
1047                 cz= (char *)&zsamp;
1048                 cz[ACOMP]= ct[0];
1049                 cz[BCOMP]= ct[1];
1050                 cz[GCOMP]= ct[2];
1051         }
1052         else if(*ct==2) {
1053                 ct= ((char *)rz);
1054                 ct+= 4+2*16*(ys & 15)+2*(xs & 15);
1055                 zsamp= *rz;
1056         
1057                 cz= (char *)&zsamp;
1058                 cz[BCOMP]= ct[0];
1059                 cz[GCOMP]= ct[1];
1060         }
1061         else if(*ct==1) {
1062                 ct= ((char *)rz);
1063                 ct+= 4+16*(ys & 15)+(xs & 15);
1064                 zsamp= *rz;
1065
1066                 cz= (char *)&zsamp;
1067                 cz[GCOMP]= ct[0];
1068
1069         }
1070         else {
1071                 /* got warning on this for 64 bits.... */
1072                 /* but it's working code! in this case rz is not a pointer but zvalue (ton) */
1073                 zsamp= GET_INT_FROM_POINTER(rz);
1074         }
1075
1076         /* tricky stuff here; we use ints which can overflow easily with bias values */
1077         
1078         if(zsamp > zs) return 1.0;              /* absolute no shadow */
1079         else if(zs < - 0x7FFFFE00 + bias) return 1.0;   /* extreme close to clipstart */
1080         else if(zsamp < zs-bias) return 0.0 ;   /* absolute in shadow */
1081         else {                                  /* soft area */
1082                 
1083                 temp=  ( (float)(zs- zsamp) )/(float)bias;
1084                 return 1.0 - temp*temp;
1085                         
1086         }
1087 }
1088
1089 /* the externally called shadow testing (reading) function */
1090 /* return 1.0: no shadow at all */
1091 float testshadowbuf(Render *re, ShadBuf *shb, float *rco, float *dxco, float *dyco, float inp, float mat_bias)
1092 {
1093         ShadSampleBuf *shsample;
1094         float fac, co[4], dx[3], dy[3], shadfac=0.0f;
1095         float xs1,ys1, siz, *jit, *weight, xres, yres, biasf;
1096         int xs, ys, zs, bias, *rz;
1097         short a, num;
1098         
1099         /* crash preventer */
1100         if(shb->buffers.first==NULL)
1101                 return 1.0f;
1102         
1103         if(inp <= 0.0f) return 0.0f;
1104
1105         /* rotate renderco en osaco */
1106         siz= 0.5f*(float)shb->size;
1107         VECCOPY(co, rco);
1108         co[3]= 1.0f;
1109
1110         Mat4MulVec4fl(shb->persmat, co);        /* rational hom co */
1111
1112         xs1= siz*(1.0f+co[0]/co[3]);
1113         ys1= siz*(1.0f+co[1]/co[3]);
1114
1115         /* Clip for z: clipsta and clipend clip values of the shadow buffer. We
1116                 * can test for -1.0/1.0 because of the properties of the
1117                 * coordinate transformations. */
1118         fac= (co[2]/co[3]);
1119
1120         if(fac>=1.0f) {
1121                 return 0.0f;
1122         } else if(fac<= -1.0f) {
1123                 return 1.0f;
1124         }
1125
1126         zs= ((float)0x7FFFFFFF)*fac;
1127
1128         /* take num*num samples, increase area with fac */
1129         num= get_render_shadow_samples(&re->r, shb->samp);
1130         num= num*num;
1131         fac= shb->soft;
1132         
1133         if(mat_bias!=0.0f) biasf= shb->bias*mat_bias;
1134         else biasf= shb->bias;
1135         /* with inp==1.0, bias is half the size. correction value was 1.1, giving errors 
1136            on cube edges, with one side being almost frontal lighted (ton)  */
1137         bias= (1.5f-inp*inp)*biasf;
1138         
1139         if(num==1) {
1140                 for(shsample= shb->buffers.first; shsample; shsample= shsample->next)
1141                         shadfac += readshadowbuf(shb, shsample, bias, (int)xs1, (int)ys1, zs);
1142                 
1143                 return shadfac/(float)shb->totbuf;
1144         }
1145
1146         /* calculate filter size */
1147         co[0]= rco[0]+dxco[0];
1148         co[1]= rco[1]+dxco[1];
1149         co[2]= rco[2]+dxco[2];
1150         co[3]= 1.0;
1151         Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
1152         dx[0]= xs1- siz*(1.0+co[0]/co[3]);
1153         dx[1]= ys1- siz*(1.0+co[1]/co[3]);
1154         
1155         co[0]= rco[0]+dyco[0];
1156         co[1]= rco[1]+dyco[1];
1157         co[2]= rco[2]+dyco[2];
1158         co[3]= 1.0;
1159         Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
1160         dy[0]= xs1- siz*(1.0+co[0]/co[3]);
1161         dy[1]= ys1- siz*(1.0+co[1]/co[3]);
1162         
1163         xres= fac*( fabs(dx[0])+fabs(dy[0]) );
1164         yres= fac*( fabs(dx[1])+fabs(dy[1]) );
1165         if(xres<fac) xres= fac;
1166         if(yres<fac) yres= fac;
1167         
1168         xs1-= (xres)/2;
1169         ys1-= (yres)/2;
1170
1171         if(xres<16.0f && yres<16.0f) {
1172                 shsample= shb->buffers.first;
1173             if(firstreadshadbuf(shb, shsample, &rz, (int)xs1, (int)ys1, 0)) {
1174                         if(firstreadshadbuf(shb, shsample, &rz, (int)(xs1+xres), (int)ys1, 1)) {
1175                                 if(firstreadshadbuf(shb, shsample, &rz, (int)xs1, (int)(ys1+yres), 1)) {
1176                                         if(firstreadshadbuf(shb, shsample, &rz, (int)(xs1+xres), (int)(ys1+yres), 1)) {
1177                                                 return readshadowbuf(shb, shsample, bias,(int)xs1, (int)ys1, zs);
1178                                         }
1179                                 }
1180                         }
1181             }
1182         }
1183         
1184         for(shsample= shb->buffers.first; shsample; shsample= shsample->next) {
1185                 jit= shb->jit;
1186                 weight= shb->weight;
1187                 
1188                 for(a=num; a>0; a--, jit+=2, weight++) {
1189                         /* instead of jit i tried random: ugly! */
1190                         /* note: the plus 0.5 gives best sampling results, jit goes from -0.5 to 0.5 */
1191                         /* xs1 and ys1 are already corrected to be corner of sample area */
1192                         xs= xs1 + xres*(jit[0] + 0.5f);
1193                         ys= ys1 + yres*(jit[1] + 0.5f);
1194                         
1195                         shadfac+= *weight * readshadowbuf(shb, shsample, bias, xs, ys, zs);
1196                 }
1197         }
1198
1199         /* Renormalizes for the sample number: */
1200         return shadfac/(float)shb->totbuf;
1201 }
1202
1203 /* different function... sampling behind clipend can be LIGHT, bias is negative! */
1204 /* return: light */
1205 static float readshadowbuf_halo(ShadBuf *shb, ShadSampleBuf *shsample, int xs, int ys, int zs)
1206 {
1207         float temp;
1208         int *rz, ofs;
1209         int bias, zbias, zsamp;
1210         char *ct, *cz;
1211
1212         /* negative! The other side is more important */
1213         bias= -shb->bias;
1214         
1215         /* simpleclip */
1216         if(xs<0 || ys<0) return 0.0;
1217         if(xs>=shb->size || ys>=shb->size) return 0.0;
1218
1219         /* calc z */
1220         ofs= (ys>>4)*(shb->size>>4) + (xs>>4);
1221         ct= shsample->cbuf+ofs;
1222         rz= *( (int **)(shsample->zbuf+ofs) );
1223
1224         if(*ct==3) {
1225                 ct= ((char *)rz)+3*16*(ys & 15)+3*(xs & 15);
1226                 cz= (char *)&zsamp;
1227                 zsamp= 0;
1228                 cz[ACOMP]= ct[0];
1229                 cz[BCOMP]= ct[1];
1230                 cz[GCOMP]= ct[2];
1231         }
1232         else if(*ct==2) {
1233                 ct= ((char *)rz);
1234                 ct+= 4+2*16*(ys & 15)+2*(xs & 15);
1235                 zsamp= *rz;
1236         
1237                 cz= (char *)&zsamp;
1238                 cz[BCOMP]= ct[0];
1239                 cz[GCOMP]= ct[1];
1240         }
1241         else if(*ct==1) {
1242                 ct= ((char *)rz);
1243                 ct+= 4+16*(ys & 15)+(xs & 15);
1244                 zsamp= *rz;
1245
1246                 cz= (char *)&zsamp;
1247                 cz[GCOMP]= ct[0];
1248
1249         }
1250         else {
1251                 /* same as before */
1252                 /* still working code! (ton) */
1253                 zsamp= GET_INT_FROM_POINTER(rz);
1254         }
1255
1256         /* NO schadow when sampled at 'eternal' distance */
1257
1258         if(zsamp >= 0x7FFFFE00) return 1.0; 
1259
1260         if(zsamp > zs) return 1.0;              /* absolute no shadww */
1261         else {
1262                 /* bias is negative, so the (zs-bias) can be beyond 0x7fffffff */
1263                 zbias= 0x7fffffff - zs;
1264                 if(zbias > -bias) {
1265                         if( zsamp < zs-bias) return 0.0 ;       /* absolute in shadow */
1266                 }
1267                 else return 0.0 ;       /* absolute shadow */
1268         }
1269
1270         /* soft area */
1271         
1272         temp=  ( (float)(zs- zsamp) )/(float)bias;
1273         return 1.0 - temp*temp;
1274 }
1275
1276
1277 float shadow_halo(LampRen *lar, float *p1, float *p2)
1278 {
1279         /* p1 p2 already are rotated in spot-space */
1280         ShadBuf *shb= lar->shb;
1281         ShadSampleBuf *shsample;
1282         float co[4], siz;
1283         float labda, labdao, labdax, labday, ldx, ldy;
1284         float zf, xf1, yf1, zf1, xf2, yf2, zf2;
1285         float count, lightcount;
1286         int x, y, z, xs1, ys1;
1287         int dx = 0, dy = 0;
1288         
1289         siz= 0.5*(float)shb->size;
1290         
1291         co[0]= p1[0];
1292         co[1]= p1[1];
1293         co[2]= p1[2]/lar->sh_zfac;
1294         co[3]= 1.0;
1295         Mat4MulVec4fl(shb->winmat, co); /* rational hom co */
1296         xf1= siz*(1.0+co[0]/co[3]);
1297         yf1= siz*(1.0+co[1]/co[3]);
1298         zf1= (co[2]/co[3]);
1299
1300
1301         co[0]= p2[0];
1302         co[1]= p2[1];
1303         co[2]= p2[2]/lar->sh_zfac;
1304         co[3]= 1.0;
1305         Mat4MulVec4fl(shb->winmat, co); /* rational hom co */
1306         xf2= siz*(1.0+co[0]/co[3]);
1307         yf2= siz*(1.0+co[1]/co[3]);
1308         zf2= (co[2]/co[3]);
1309
1310         /* the 2dda (a pixel line formula) */
1311
1312         xs1= (int)xf1;
1313         ys1= (int)yf1;
1314
1315         if(xf1 != xf2) {
1316                 if(xf2-xf1 > 0.0) {
1317                         labdax= (xf1-xs1-1.0)/(xf1-xf2);
1318                         ldx= -shb->shadhalostep/(xf1-xf2);
1319                         dx= shb->shadhalostep;
1320                 }
1321                 else {
1322                         labdax= (xf1-xs1)/(xf1-xf2);
1323                         ldx= shb->shadhalostep/(xf1-xf2);
1324                         dx= -shb->shadhalostep;
1325                 }
1326         }
1327         else {
1328                 labdax= 1.0;
1329                 ldx= 0.0;
1330         }
1331
1332         if(yf1 != yf2) {
1333                 if(yf2-yf1 > 0.0) {
1334                         labday= (yf1-ys1-1.0)/(yf1-yf2);
1335                         ldy= -shb->shadhalostep/(yf1-yf2);
1336                         dy= shb->shadhalostep;
1337                 }
1338                 else {
1339                         labday= (yf1-ys1)/(yf1-yf2);
1340                         ldy= shb->shadhalostep/(yf1-yf2);
1341                         dy= -shb->shadhalostep;
1342                 }
1343         }
1344         else {
1345                 labday= 1.0;
1346                 ldy= 0.0;
1347         }
1348         
1349         x= xs1;
1350         y= ys1;
1351         labda= count= lightcount= 0.0;
1352
1353 /* printf("start %x %x  \n", (int)(0x7FFFFFFF*zf1), (int)(0x7FFFFFFF*zf2)); */
1354
1355         while(1) {
1356                 labdao= labda;
1357                 
1358                 if(labdax==labday) {
1359                         labdax+= ldx;
1360                         x+= dx;
1361                         labday+= ldy;
1362                         y+= dy;
1363                 }
1364                 else {
1365                         if(labdax<labday) {
1366                                 labdax+= ldx;
1367                                 x+= dx;
1368                         } else {
1369                                 labday+= ldy;
1370                                 y+= dy;
1371                         }
1372                 }
1373                 
1374                 labda= MIN2(labdax, labday);
1375                 if(labda==labdao || labda>=1.0) break;
1376                 
1377                 zf= zf1 + labda*(zf2-zf1);
1378                 count+= (float)shb->totbuf;
1379
1380                 if(zf<= -1.0) lightcount += 1.0;        /* close to the spot */
1381                 else {
1382                 
1383                         /* make sure, behind the clipend we extend halolines. */
1384                         if(zf>=1.0) z= 0x7FFFF000;
1385                         else z= (int)(0x7FFFF000*zf);
1386                         
1387                         for(shsample= shb->buffers.first; shsample; shsample= shsample->next)
1388                                 lightcount+= readshadowbuf_halo(shb, shsample, x, y, z);
1389                         
1390                 }
1391         }
1392         
1393         if(count!=0.0) return (lightcount/count);
1394         return 0.0;
1395         
1396 }
1397
1398
1399 /* ********************* Irregular Shadow Buffer (ISB) ************* */
1400 /* ********** storage of all view samples in a raster of lists ***** */
1401
1402 /* based on several articles describing this method, like:
1403 The Irregular Z-Buffer and its Application to Shadow Mapping
1404 Gregory S. Johnson - William R. Mark - Christopher A. Burns 
1405 and
1406 Alias-Free Shadow Maps
1407 Timo Aila and Samuli Laine
1408 */
1409
1410 /* bsp structure (actually kd tree) */
1411
1412 #define BSPMAX_SAMPLE   128
1413 #define BSPMAX_DEPTH    32
1414
1415 /* aligned with struct rctf */
1416 typedef struct Boxf {
1417         float xmin, xmax;
1418         float ymin, ymax;
1419         float zmin, zmax;
1420 } Boxf;
1421
1422 typedef struct ISBBranch {
1423         struct ISBBranch *left, *right;
1424         float divider[2];
1425         Boxf box;
1426         short totsamp, index, full, unused;
1427         ISBSample **samples;
1428 } ISBBranch;
1429
1430 typedef struct BSPFace {
1431         Boxf box;
1432         float *v1, *v2, *v3, *v4;
1433         int obi;                /* object for face lookup */
1434         int facenr;             /* index to retrieve VlakRen */
1435         int type;               /* only for strand now */
1436         short shad_alpha, is_full;
1437         
1438         /* strand caching data, optimize for point_behind_strand() */
1439         float radline, radline_end, len;
1440         float vec1[3], vec2[3], rc[3];
1441 } BSPFace;
1442
1443 /* boxes are in lamp projection */
1444 static void init_box(Boxf *box)
1445 {
1446         box->xmin= 1000000.0f;
1447         box->xmax= 0;
1448         box->ymin= 1000000.0f;
1449         box->ymax= 0;
1450         box->zmin= 0x7FFFFFFF;
1451         box->zmax= - 0x7FFFFFFF;
1452 }
1453
1454 /* use v1 to calculate boundbox */
1455 static void bound_boxf(Boxf *box, float *v1)
1456 {
1457         if(v1[0] < box->xmin) box->xmin= v1[0];
1458         if(v1[0] > box->xmax) box->xmax= v1[0];
1459         if(v1[1] < box->ymin) box->ymin= v1[1];
1460         if(v1[1] > box->ymax) box->ymax= v1[1];
1461         if(v1[2] < box->zmin) box->zmin= v1[2];
1462         if(v1[2] > box->zmax) box->zmax= v1[2];
1463 }
1464
1465 /* use v1 to calculate boundbox */
1466 static void bound_rectf(rctf *box, float *v1)
1467 {
1468         if(v1[0] < box->xmin) box->xmin= v1[0];
1469         if(v1[0] > box->xmax) box->xmax= v1[0];
1470         if(v1[1] < box->ymin) box->ymin= v1[1];
1471         if(v1[1] > box->ymax) box->ymax= v1[1];
1472 }
1473
1474
1475 /* halfway splitting, for initializing a more regular tree */
1476 static void isb_bsp_split_init(ISBBranch *root, MemArena *mem, int level)
1477 {
1478         
1479         /* if level > 0 we create new branches and go deeper*/
1480         if(level > 0) {
1481                 ISBBranch *left, *right;
1482                 int i;
1483                 
1484                 /* splitpoint */
1485                 root->divider[0]= 0.5f*(root->box.xmin+root->box.xmax);
1486                 root->divider[1]= 0.5f*(root->box.ymin+root->box.ymax);
1487                 
1488                 /* find best splitpoint */
1489                 if(root->box.xmax-root->box.xmin > root->box.ymax-root->box.ymin)
1490                         i= root->index= 0;
1491                 else
1492                         i= root->index= 1;
1493                 
1494                 left= root->left= BLI_memarena_alloc(mem, sizeof(ISBBranch));
1495                 right= root->right= BLI_memarena_alloc(mem, sizeof(ISBBranch));
1496                 
1497                 /* box info */
1498                 left->box= root->box;
1499                 right->box= root->box;
1500                 if(i==0) {
1501                         left->box.xmax= root->divider[0];
1502                         right->box.xmin= root->divider[0];
1503                 }
1504                 else {
1505                         left->box.ymax= root->divider[1];
1506                         right->box.ymin= root->divider[1];
1507                 }
1508                 isb_bsp_split_init(left, mem, level-1);
1509                 isb_bsp_split_init(right, mem, level-1);
1510         }
1511         else {
1512                 /* we add sample array */
1513                 root->samples= BLI_memarena_alloc(mem, BSPMAX_SAMPLE*sizeof(void *));
1514         }
1515 }
1516
1517 /* note; if all samples on same location we just spread them over 2 new branches */
1518 static void isb_bsp_split(ISBBranch *root, MemArena *mem)
1519 {
1520         ISBBranch *left, *right;
1521         ISBSample *samples[BSPMAX_SAMPLE];
1522         int a, i;
1523
1524         /* splitpoint */
1525         root->divider[0]= root->divider[1]= 0.0f;
1526         for(a=BSPMAX_SAMPLE-1; a>=0; a--) {
1527                 root->divider[0]+= root->samples[a]->zco[0];
1528                 root->divider[1]+= root->samples[a]->zco[1];
1529         }
1530         root->divider[0]/= BSPMAX_SAMPLE;
1531         root->divider[1]/= BSPMAX_SAMPLE;
1532         
1533         /* find best splitpoint */
1534         if(root->box.xmax-root->box.xmin > root->box.ymax-root->box.ymin)
1535                 i= root->index= 0;
1536         else
1537                 i= root->index= 1;
1538         
1539         /* new branches */
1540         left= root->left= BLI_memarena_alloc(mem, sizeof(ISBBranch));
1541         right= root->right= BLI_memarena_alloc(mem, sizeof(ISBBranch));
1542
1543         /* new sample array */
1544         left->samples= BLI_memarena_alloc(mem, BSPMAX_SAMPLE*sizeof(void *));
1545         right->samples= samples; // tmp
1546         
1547         /* split samples */
1548         for(a=BSPMAX_SAMPLE-1; a>=0; a--) {
1549                 int comp= 0;
1550                 /* this prevents adding samples all to 1 branch when divider is equal to samples */
1551                 if(root->samples[a]->zco[i] == root->divider[i])
1552                         comp= a & 1;
1553                 else if(root->samples[a]->zco[i] < root->divider[i])
1554                         comp= 1;
1555                 
1556                 if(comp==1) {
1557                         left->samples[left->totsamp]= root->samples[a];
1558                         left->totsamp++;
1559                 }
1560                 else {
1561                         right->samples[right->totsamp]= root->samples[a];
1562                         right->totsamp++;
1563                 }
1564         }
1565         
1566         /* copy samples from tmp */
1567         memcpy(root->samples, samples, right->totsamp*(sizeof(void *)));
1568         right->samples= root->samples;
1569         root->samples= NULL;
1570         
1571         /* box info */
1572         left->box= root->box;
1573         right->box= root->box;
1574         if(i==0) {
1575                 left->box.xmax= root->divider[0];
1576                 right->box.xmin= root->divider[0];
1577         }
1578         else {
1579                 left->box.ymax= root->divider[1];
1580                 right->box.ymin= root->divider[1];
1581         }
1582 }
1583
1584 /* inserts sample in main tree, also splits on threshold */
1585 /* returns 1 if error */
1586 static int isb_bsp_insert(ISBBranch *root, MemArena *memarena, ISBSample *sample)
1587 {
1588         ISBBranch *bspn= root;
1589         float *zco= sample->zco;
1590         int i= 0;
1591         
1592         /* debug counter, also used to check if something was filled in ever */
1593         root->totsamp++;
1594         
1595         /* going over branches until last one found */
1596         while(bspn->left) {
1597                 if(zco[bspn->index] <= bspn->divider[bspn->index])
1598                         bspn= bspn->left;
1599                 else
1600                         bspn= bspn->right;
1601                 i++;
1602         }
1603         /* bspn now is the last branch */
1604         
1605         if(bspn->totsamp==BSPMAX_SAMPLE) {
1606                 printf("error in bsp branch\n");        /* only for debug, cannot happen */
1607                 return 1;
1608         }
1609         
1610         /* insert */
1611         bspn->samples[bspn->totsamp]= sample;
1612         bspn->totsamp++;
1613
1614         /* split if allowed and needed */
1615         if(bspn->totsamp==BSPMAX_SAMPLE) {
1616                 if(i==BSPMAX_DEPTH) {
1617                         bspn->totsamp--;        /* stop filling in... will give errors */
1618                         return 1;
1619                 }
1620                 isb_bsp_split(bspn, memarena);
1621         }
1622         return 0;
1623 }
1624
1625 static float VecLen2f( float *v1, float *v2)
1626 {
1627         float x= v1[0]-v2[0];
1628         float y= v1[1]-v2[1];
1629         return (float)sqrt(x*x+y*y);
1630 }
1631
1632 /* initialize vars in face, for optimal point-in-face test */
1633 static void bspface_init_strand(BSPFace *face) 
1634 {
1635         
1636         face->radline= 0.5f*VecLen2f(face->v1, face->v2);
1637         
1638         VecMidf(face->vec1, face->v1, face->v2);
1639         if(face->v4)
1640                 VecMidf(face->vec2, face->v3, face->v4);
1641         else
1642                 VECCOPY(face->vec2, face->v3);
1643         
1644         face->rc[0]= face->vec2[0]-face->vec1[0];
1645         face->rc[1]= face->vec2[1]-face->vec1[1];
1646         face->rc[2]= face->vec2[2]-face->vec1[2];
1647         
1648         face->len= face->rc[0]*face->rc[0]+ face->rc[1]*face->rc[1];
1649         
1650         if(face->len!=0.0f) {
1651                 face->radline_end= face->radline/sqrt(face->len);
1652                 face->len= 1.0f/face->len;
1653         }
1654 }
1655
1656 /* brought back to a simple 2d case */
1657 static int point_behind_strand(float *p, BSPFace *face)
1658 {
1659         /* v1 - v2 is radius, v1 - v3 length */
1660         float dist, rc[2], pt[2];
1661         
1662         /* using code from PdistVL2Dfl(), distance vec to line-piece */
1663
1664         if(face->len==0.0f) {
1665                 rc[0]= p[0]-face->vec1[0];
1666                 rc[1]= p[1]-face->vec1[1];
1667                 dist= (float)(sqrt(rc[0]*rc[0]+ rc[1]*rc[1]));
1668                 
1669                 if(dist < face->radline)
1670                         return 1;
1671         }
1672         else {
1673                 float labda= ( face->rc[0]*(p[0]-face->vec1[0]) + face->rc[1]*(p[1]-face->vec1[1]) )*face->len;
1674                 
1675                 if(labda > -face->radline_end && labda < 1.0f+face->radline_end) { 
1676                         /* hesse for dist: */
1677                         //dist= (float)(fabs( (p[0]-vec2[0])*rc[1] + (p[1]-vec2[1])*rc[0])/len);
1678                         
1679                         pt[0]= labda*face->rc[0]+face->vec1[0];
1680                         pt[1]= labda*face->rc[1]+face->vec1[1];
1681                         
1682                         rc[0]= pt[0]-p[0];
1683                         rc[1]= pt[1]-p[1];
1684                         dist= (float)sqrt(rc[0]*rc[0]+ rc[1]*rc[1]);
1685                         
1686                         if(dist < face->radline) {
1687                                 float zval= face->vec1[2] + labda*face->rc[2];
1688                                 if(p[2] > zval)
1689                                         return 1;
1690                         }
1691                 }
1692         }
1693         return 0;
1694 }
1695
1696
1697 /* return 1 if inside. code derived from src/parametrizer.c */
1698 static int point_behind_tria2d(float *p, float *v1, float *v2, float *v3)
1699 {
1700         float a[2], c[2], h[2], div;
1701         float u, v;
1702         
1703         a[0] = v2[0] - v1[0];
1704         a[1] = v2[1] - v1[1];
1705         c[0] = v3[0] - v1[0];
1706         c[1] = v3[1] - v1[1];
1707         
1708         div = a[0]*c[1] - a[1]*c[0];
1709         if(div==0.0f)
1710                 return 0;
1711         
1712         h[0] = p[0] - v1[0];
1713         h[1] = p[1] - v1[1];
1714         
1715         div = 1.0f/div;
1716         
1717         u = (h[0]*c[1] - h[1]*c[0])*div;
1718         if(u >= 0.0f) {
1719                 v = (a[0]*h[1] - a[1]*h[0])*div;
1720                 if(v >= 0.0f) {
1721                         if( u + v <= 1.0f) {
1722                                 /* inside, now check if point p is behind */
1723                                 float z=  (1.0f-u-v)*v1[2] + u*v2[2] + v*v3[2];
1724                                 if(z <= p[2])
1725                                         return 1;
1726                         }
1727                 }
1728         }
1729         
1730         return 0;
1731 }
1732
1733 #if 0
1734 /* tested these calls, but it gives inaccuracy, 'side' cannot be found reliably using v3 */
1735
1736 /* check if line v1-v2 has all rect points on other side of point v3 */
1737 static int rect_outside_line(rctf *rect, float *v1, float *v2, float *v3)
1738 {
1739         float a, b, c;
1740         int side;
1741         
1742         /* line formula for v1-v2 */
1743         a= v2[1]-v1[1];
1744         b= v1[0]-v2[0];
1745         c= -a*v1[0] - b*v1[1];
1746         side= a*v3[0] + b*v3[1] + c < 0.0f;
1747         
1748         /* the four quad points */
1749         if( side==(rect->xmin*a + rect->ymin*b + c >= 0.0f) )
1750                 if( side==(rect->xmax*a + rect->ymin*b + c >= 0.0f) )
1751                         if( side==(rect->xmax*a + rect->ymax*b + c >= 0.0f) )
1752                                 if( side==(rect->xmin*a + rect->ymax*b + c >= 0.0f) )
1753                                         return 1;
1754         return 0;
1755 }
1756
1757 /* check if one of the triangle edges separates all rect points on 1 side */
1758 static int rect_isect_tria(rctf *rect, float *v1, float *v2, float *v3)
1759 {
1760         if(rect_outside_line(rect, v1, v2, v3))
1761                 return 0;
1762         if(rect_outside_line(rect, v2, v3, v1))
1763                 return 0;
1764         if(rect_outside_line(rect, v3, v1, v2))
1765                 return 0;
1766         return 1;
1767 }
1768 #endif
1769
1770 /* if face overlaps a branch, it executes func. recursive */
1771 static void isb_bsp_face_inside(ISBBranch *bspn, BSPFace *face)
1772 {
1773         
1774         /* are we descending? */
1775         if(bspn->left) {
1776                 /* hrmf, the box struct cannot be addressed with index */
1777                 if(bspn->index==0) {
1778                         if(face->box.xmin <= bspn->divider[0])
1779                                 isb_bsp_face_inside(bspn->left, face);
1780                         if(face->box.xmax > bspn->divider[0])
1781                                 isb_bsp_face_inside(bspn->right, face);
1782                 }
1783                 else {
1784                         if(face->box.ymin <= bspn->divider[1])
1785                                 isb_bsp_face_inside(bspn->left, face);
1786                         if(face->box.ymax > bspn->divider[1])
1787                                 isb_bsp_face_inside(bspn->right, face);
1788                 }
1789         }
1790         else {
1791                 /* else: end branch reached */
1792                 int a;
1793                 
1794                 if(bspn->totsamp==0) return;
1795                 
1796                 /* check for nodes entirely in shadow, can be skipped */
1797                 if(bspn->totsamp==bspn->full)
1798                         return;
1799                 
1800                 /* if bsp node is entirely in front of face, give up */
1801                 if(bspn->box.zmax < face->box.zmin)
1802                         return;
1803                 
1804                 /* if face boundbox is outside of branch rect, give up */
1805                 if(0==BLI_isect_rctf((rctf *)&face->box, (rctf *)&bspn->box, NULL))
1806                         return;
1807                 
1808                 /* test all points inside branch */
1809                 for(a=bspn->totsamp-1; a>=0; a--) {
1810                         ISBSample *samp= bspn->samples[a];
1811                         
1812                         if((samp->facenr!=face->facenr || samp->obi!=face->obi) && samp->shadfac) {
1813                                 if(face->box.zmin < samp->zco[2]) {
1814                                         if(BLI_in_rctf((rctf *)&face->box, samp->zco[0], samp->zco[1])) {
1815                                                 int inshadow= 0;
1816                                                 
1817                                                 if(face->type) {
1818                                                         if(point_behind_strand(samp->zco, face)) 
1819                                                                 inshadow= 1;
1820                                                 }
1821                                                 else if( point_behind_tria2d(samp->zco, face->v1, face->v2, face->v3))
1822                                                         inshadow= 1;
1823                                                 else if(face->v4 && point_behind_tria2d(samp->zco, face->v1, face->v3, face->v4))
1824                                                         inshadow= 1;
1825
1826                                                 if(inshadow) {
1827                                                         *(samp->shadfac) += face->shad_alpha;
1828                                                         /* optimize; is_full means shad_alpha==4096 */
1829                                                         if(*(samp->shadfac) >= 4096 || face->is_full) {
1830                                                                 bspn->full++;
1831                                                                 samp->shadfac= NULL;
1832                                                         }
1833                                                 }
1834                                         }
1835                                 }
1836                         }
1837                 }
1838         }
1839 }
1840
1841 /* based on available samples, recalculate the bounding box for bsp nodes, recursive */
1842 static void isb_bsp_recalc_box(ISBBranch *root)
1843 {
1844         if(root->left) {
1845                 isb_bsp_recalc_box(root->left);
1846                 isb_bsp_recalc_box(root->right);
1847         }
1848         else if(root->totsamp) {
1849                 int a;
1850                 
1851                 init_box(&root->box);
1852                 for(a=root->totsamp-1; a>=0; a--)
1853                         bound_boxf(&root->box, root->samples[a]->zco);
1854         }       
1855 }
1856
1857 /* callback function for zbuf clip */
1858 static void isb_bsp_test_strand(ZSpan *zspan, int obi, int zvlnr, float *v1, float *v2, float *v3, float *v4)
1859 {
1860         BSPFace face;
1861         
1862         face.v1= v1;
1863         face.v2= v2;
1864         face.v3= v3;
1865         face.v4= v4;
1866         face.obi= obi;
1867         face.facenr= zvlnr & ~RE_QUAD_OFFS;
1868         face.type= R_STRAND;
1869         if(R.osa)
1870                 face.shad_alpha= (short)ceil(4096.0f*zspan->shad_alpha/(float)R.osa);
1871         else
1872                 face.shad_alpha= (short)ceil(4096.0f*zspan->shad_alpha);
1873         
1874         face.is_full= (zspan->shad_alpha==1.0f);
1875         
1876         /* setup boundbox */
1877         init_box(&face.box);
1878         bound_boxf(&face.box, v1);
1879         bound_boxf(&face.box, v2);
1880         bound_boxf(&face.box, v3);
1881         if(v4)
1882                 bound_boxf(&face.box, v4);
1883         
1884         /* optimize values */
1885         bspface_init_strand(&face);
1886         
1887         isb_bsp_face_inside((ISBBranch *)zspan->rectz, &face);
1888         
1889 }
1890
1891 /* callback function for zbuf clip */
1892 static void isb_bsp_test_face(ZSpan *zspan, int obi, int zvlnr, float *v1, float *v2, float *v3, float *v4) 
1893 {
1894         BSPFace face;
1895         
1896         face.v1= v1;
1897         face.v2= v2;
1898         face.v3= v3;
1899         face.v4= v4;
1900         face.obi= obi;
1901         face.facenr= zvlnr & ~RE_QUAD_OFFS;
1902         face.type= 0;
1903         if(R.osa)
1904                 face.shad_alpha= (short)ceil(4096.0f*zspan->shad_alpha/(float)R.osa);
1905         else
1906                 face.shad_alpha= (short)ceil(4096.0f*zspan->shad_alpha);
1907         
1908         face.is_full= (zspan->shad_alpha==1.0f);
1909         
1910         /* setup boundbox */
1911         init_box(&face.box);
1912         bound_boxf(&face.box, v1);
1913         bound_boxf(&face.box, v2);
1914         bound_boxf(&face.box, v3);
1915         if(v4)
1916                 bound_boxf(&face.box, v4);
1917
1918         isb_bsp_face_inside((ISBBranch *)zspan->rectz, &face);
1919 }
1920
1921 static int testclip_minmax(float *ho, float *minmax)
1922 {
1923         float wco= ho[3];
1924         int flag= 0;
1925         
1926         if( ho[0] > minmax[1]*wco) flag = 1;
1927         else if( ho[0]< minmax[0]*wco) flag = 2;
1928         
1929         if( ho[1] > minmax[3]*wco) flag |= 4;
1930         else if( ho[1]< minmax[2]*wco) flag |= 8;
1931         
1932         return flag;
1933 }
1934
1935 /* main loop going over all faces and check in bsp overlaps, fill in shadfac values */
1936 static void isb_bsp_fillfaces(Render *re, LampRen *lar, ISBBranch *root)
1937 {
1938         ObjectInstanceRen *obi;
1939         ObjectRen *obr;
1940         ShadBuf *shb= lar->shb;
1941         ZSpan zspan, zspanstrand;
1942         VlakRen *vlr= NULL;
1943         Material *ma= NULL;
1944         float minmaxf[4], winmat[4][4];
1945         int size= shb->size;
1946         int i, a, ok=1, lay= -1;
1947         
1948         /* further optimize, also sets minz maxz */
1949         isb_bsp_recalc_box(root);
1950         
1951         /* extra clipping for minmax */
1952         minmaxf[0]= (2.0f*root->box.xmin - size-2.0f)/size;
1953         minmaxf[1]= (2.0f*root->box.xmax - size+2.0f)/size;
1954         minmaxf[2]= (2.0f*root->box.ymin - size-2.0f)/size;
1955         minmaxf[3]= (2.0f*root->box.ymax - size+2.0f)/size;
1956         
1957         if(lar->mode & (LA_LAYER|LA_LAYER_SHADOW)) lay= lar->lay;
1958         
1959         /* (ab)use zspan, since we use zbuffer clipping code */
1960         zbuf_alloc_span(&zspan, size, size, re->clipcrop);
1961         
1962         zspan.zmulx=  ((float)size)/2.0f;
1963         zspan.zmuly=  ((float)size)/2.0f;
1964         zspan.zofsx= -0.5f;
1965         zspan.zofsy= -0.5f;
1966         
1967         /* pass on bsp root to zspan */
1968         zspan.rectz= (int *)root;
1969         
1970         /* filling methods */
1971         zspanstrand= zspan;
1972         //      zspan.zbuflinefunc= zbufline_onlyZ;
1973         zspan.zbuffunc= isb_bsp_test_face;
1974         zspanstrand.zbuffunc= isb_bsp_test_strand;
1975         
1976         for(i=0, obi=re->instancetable.first; obi; i++, obi=obi->next) {
1977                 obr= obi->obr;
1978
1979                 if(obi->flag & R_TRANSFORMED)
1980                         Mat4MulMat4(winmat, obi->mat, shb->persmat);
1981                 else
1982                         Mat4CpyMat4(winmat, shb->persmat);
1983
1984                 for(a=0; a<obr->totvlak; a++) {
1985                         
1986                         if((a & 255)==0) vlr= obr->vlaknodes[a>>8].vlak;
1987                         else vlr++;
1988                         
1989                         /* note, these conditions are copied in shadowbuf_autoclip() */
1990                         if(vlr->mat!= ma) {
1991                                 ma= vlr->mat;
1992                                 ok= 1;
1993                                 if((ma->mode & MA_SHADBUF)==0) ok= 0;
1994                                 if(ma->material_type == MA_TYPE_WIRE) ok= 0;
1995                                 zspanstrand.shad_alpha= zspan.shad_alpha= ma->shad_alpha;
1996                         }
1997                         
1998                         if(ok && (obi->lay & lay)) {
1999                                 float hoco[4][4];
2000                                 int c1, c2, c3, c4=0;
2001                                 int d1, d2, d3, d4=0;
2002                                 int partclip;
2003                                 
2004                                 /* create hocos per face, it is while render */
2005                                 projectvert(vlr->v1->co, winmat, hoco[0]); d1= testclip_minmax(hoco[0], minmaxf);
2006                                 projectvert(vlr->v2->co, winmat, hoco[1]); d2= testclip_minmax(hoco[1], minmaxf);
2007                                 projectvert(vlr->v3->co, winmat, hoco[2]); d3= testclip_minmax(hoco[2], minmaxf);
2008                                 if(vlr->v4) {
2009                                         projectvert(vlr->v4->co, winmat, hoco[3]); d4= testclip_minmax(hoco[3], minmaxf);
2010                                 }
2011
2012                                 /* minmax clipping */
2013                                 if(vlr->v4) partclip= d1 & d2 & d3 & d4;
2014                                 else partclip= d1 & d2 & d3;
2015                                 
2016                                 if(partclip==0) {
2017                                         
2018                                         /* window clipping */
2019                                         c1= testclip(hoco[0]); 
2020                                         c2= testclip(hoco[1]); 
2021                                         c3= testclip(hoco[2]); 
2022                                         if(vlr->v4)
2023                                                 c4= testclip(hoco[3]); 
2024                                         
2025                                         /* ***** NO WIRE YET */                 
2026                                         if(ma->material_type == MA_TYPE_WIRE) {
2027                                                 if(vlr->v4)
2028                                                         zbufclipwire(&zspan, i, a+1, vlr->ec, hoco[0], hoco[1], hoco[2], hoco[3], c1, c2, c3, c4);
2029                                                 else
2030                                                         zbufclipwire(&zspan, i, a+1, vlr->ec, hoco[0], hoco[1], hoco[2], 0, c1, c2, c3, 0);
2031                                         }
2032                                         else if(vlr->v4) {
2033                                                 if(vlr->flag & R_STRAND)
2034                                                         zbufclip4(&zspanstrand, i, a+1, hoco[0], hoco[1], hoco[2], hoco[3], c1, c2, c3, c4);
2035                                                 else
2036                                                         zbufclip4(&zspan, i, a+1, hoco[0], hoco[1], hoco[2], hoco[3], c1, c2, c3, c4);
2037                                         }
2038                                         else
2039                                                 zbufclip(&zspan, i, a+1, hoco[0], hoco[1], hoco[2], c1, c2, c3);
2040                                         
2041                                 }
2042                         }
2043                 }
2044         }
2045         
2046         zbuf_free_span(&zspan);
2047 }
2048
2049 /* returns 1 when the viewpixel is visible in lampbuffer */
2050 static int viewpixel_to_lampbuf(ShadBuf *shb, ObjectInstanceRen *obi, VlakRen *vlr, float x, float y, float *co)
2051 {
2052         float hoco[4], v1[3], nor[3];
2053         float dface, fac, siz;
2054         
2055         RE_vlakren_get_normal(&R, obi, vlr, nor);
2056         VECCOPY(v1, vlr->v1->co);
2057         if(obi->flag & R_TRANSFORMED)
2058                 Mat4MulVecfl(obi->mat, v1);
2059
2060         /* from shadepixel() */
2061         dface= v1[0]*nor[0] + v1[1]*nor[1] + v1[2]*nor[2];
2062         hoco[3]= 1.0f;
2063         
2064         /* ortho viewplane cannot intersect using view vector originating in (0,0,0) */
2065         if(R.r.mode & R_ORTHO) {
2066                 /* x and y 3d coordinate can be derived from pixel coord and winmat */
2067                 float fx= 2.0/(R.winx*R.winmat[0][0]);
2068                 float fy= 2.0/(R.winy*R.winmat[1][1]);
2069                 
2070                 hoco[0]= (x - 0.5*R.winx)*fx - R.winmat[3][0]/R.winmat[0][0];
2071                 hoco[1]= (y - 0.5*R.winy)*fy - R.winmat[3][1]/R.winmat[1][1];
2072                 
2073                 /* using a*x + b*y + c*z = d equation, (a b c) is normal */
2074                 if(nor[2]!=0.0f)
2075                         hoco[2]= (dface - nor[0]*hoco[0] - nor[1]*hoco[1])/nor[2];
2076                 else
2077                         hoco[2]= 0.0f;
2078         }
2079         else {
2080                 float div, view[3];
2081                 
2082                 calc_view_vector(view, x, y);
2083                 
2084                 div= nor[0]*view[0] + nor[1]*view[1] + nor[2]*view[2];
2085                 if (div==0.0f) 
2086                         return 0;
2087                 
2088                 fac= dface/div;
2089                 
2090                 hoco[0]= fac*view[0];
2091                 hoco[1]= fac*view[1];
2092                 hoco[2]= fac*view[2];
2093         }
2094         
2095         /* move 3d vector to lampbuf */
2096         Mat4MulVec4fl(shb->persmat, hoco);      /* rational hom co */
2097         
2098         /* clip We can test for -1.0/1.0 because of the properties of the
2099          * coordinate transformations. */
2100         fac= fabs(hoco[3]);
2101         if(hoco[0]<-fac || hoco[0]>fac)
2102                 return 0;
2103         if(hoco[1]<-fac || hoco[1]>fac)
2104                 return 0;
2105         if(hoco[2]<-fac || hoco[2]>fac)
2106                 return 0;
2107         
2108         siz= 0.5f*(float)shb->size;
2109         co[0]= siz*(1.0f+hoco[0]/hoco[3]) -0.5f;
2110         co[1]= siz*(1.0f+hoco[1]/hoco[3]) -0.5f;
2111         co[2]= ((float)0x7FFFFFFF)*(hoco[2]/hoco[3]);
2112         
2113         /* XXXX bias, much less than normal shadbuf, or do we need a constant? */
2114         co[2] -= 0.05f*shb->bias;
2115         
2116         return 1;
2117 }
2118
2119 /* storage of shadow results, solid osa and transp case */
2120 static void isb_add_shadfac(ISBShadfacA **isbsapp, MemArena *mem, int obi, int facenr, short shadfac, short samples)
2121 {
2122         ISBShadfacA *new;
2123         float shadfacf;
2124         
2125         /* in osa case, the samples were filled in with factor 1.0/R.osa. if fewer samples we have to correct */
2126         if(R.osa)
2127                 shadfacf= ((float)shadfac*R.osa)/(4096.0*samples);
2128         else
2129                 shadfacf= ((float)shadfac)/(4096.0);
2130         
2131         new= BLI_memarena_alloc(mem, sizeof(ISBShadfacA));
2132         new->obi= obi;
2133         new->facenr= facenr & ~RE_QUAD_OFFS;
2134         new->shadfac= shadfacf;
2135         if(*isbsapp)
2136                 new->next= (*isbsapp);
2137         else
2138                 new->next= NULL;
2139         
2140         *isbsapp= new;
2141 }
2142
2143 /* adding samples, solid case */
2144 static int isb_add_samples(RenderPart *pa, ISBBranch *root, MemArena *memarena, ISBSample **samplebuf)
2145 {
2146         int xi, yi, *xcos, *ycos;
2147         int sample, bsp_err= 0;
2148         
2149         /* bsp split doesn't like to handle regular sequenes */
2150         xcos= MEM_mallocN( pa->rectx*sizeof(int), "xcos");
2151         ycos= MEM_mallocN( pa->recty*sizeof(int), "ycos");
2152         for(xi=0; xi<pa->rectx; xi++)
2153                 xcos[xi]= xi;
2154         for(yi=0; yi<pa->recty; yi++)
2155                 ycos[yi]= yi;
2156         BLI_array_randomize(xcos, sizeof(int), pa->rectx, 12345);
2157         BLI_array_randomize(ycos, sizeof(int), pa->recty, 54321);
2158         
2159         for(sample=0; sample<(R.osa?R.osa:1); sample++) {
2160                 ISBSample *samp= samplebuf[sample], *samp1;
2161                 
2162                 for(yi=0; yi<pa->recty; yi++) {
2163                         int y= ycos[yi];
2164                         for(xi=0; xi<pa->rectx; xi++) {
2165                                 int x= xcos[xi];
2166                                 samp1= samp + y*pa->rectx + x;
2167                                 if(samp1->facenr)
2168                                         bsp_err |= isb_bsp_insert(root, memarena, samp1);
2169                         }
2170                         if(bsp_err) break;
2171                 }
2172         }       
2173         
2174         MEM_freeN(xcos);
2175         MEM_freeN(ycos);
2176
2177         return bsp_err;
2178 }
2179
2180 /* solid version */
2181 /* lar->shb, pa->rectz and pa->rectp should exist */
2182 static void isb_make_buffer(RenderPart *pa, LampRen *lar)
2183 {
2184         ShadBuf *shb= lar->shb;
2185         ISBData *isbdata;
2186         ISBSample *samp, *samplebuf[16];        /* should be RE_MAX_OSA */
2187         ISBBranch root;
2188         MemArena *memarena;
2189         intptr_t *rd;
2190         int *recto, *rectp, x, y, sindex, sample, bsp_err=0;
2191         
2192         /* storage for shadow, per thread */
2193         isbdata= shb->isb_result[pa->thread];
2194         
2195         /* to map the shi->xs and ys coordinate */
2196         isbdata->minx= pa->disprect.xmin;
2197         isbdata->miny= pa->disprect.ymin;
2198         isbdata->rectx= pa->rectx;
2199         isbdata->recty= pa->recty;
2200         
2201         /* branches are added using memarena (32k branches) */
2202         memarena = BLI_memarena_new(0x8000 * sizeof(ISBBranch));
2203         BLI_memarena_use_calloc(memarena);
2204         
2205         /* samplebuf is in camera view space (pixels) */
2206         for(sample=0; sample<(R.osa?R.osa:1); sample++)
2207                 samplebuf[sample]= MEM_callocN(sizeof(ISBSample)*pa->rectx*pa->recty, "isb samplebuf");
2208         
2209         /* for end result, ISBSamples point to this in non OSA case, otherwise to pixstruct->shadfac */
2210         if(R.osa==0)
2211                 isbdata->shadfacs= MEM_callocN(pa->rectx*pa->recty*sizeof(short), "isb shadfacs");
2212         
2213         /* setup bsp root */
2214         memset(&root, 0, sizeof(ISBBranch));
2215         root.box.xmin= (float)shb->size;
2216         root.box.ymin= (float)shb->size;
2217         
2218         /* create the sample buffers */
2219         for(sindex=0, y=0; y<pa->recty; y++) {
2220                 for(x=0; x<pa->rectx; x++, sindex++) {
2221                         
2222                         /* this makes it a long function, but splitting it out would mean 10+ arguments */
2223                         /* first check OSA case */
2224                         if(R.osa) {
2225                                 rd= pa->rectdaps + sindex;
2226                                 if(*rd) {
2227                                         float xs= (float)(x + pa->disprect.xmin);
2228                                         float ys= (float)(y + pa->disprect.ymin);
2229                                         
2230                                         for(sample=0; sample<R.osa; sample++) {
2231                                                 PixStr *ps= (PixStr *)(*rd);
2232                                                 int mask= (1<<sample);
2233                                                 
2234                                                 while(ps) {
2235                                                         if(ps->mask & mask)
2236                                                                 break;
2237                                                         ps= ps->next;
2238                                                 }
2239                                                 if(ps && ps->facenr>0) {
2240                                                         ObjectInstanceRen *obi= &R.objectinstance[ps->obi];
2241                                                         ObjectRen *obr= obi->obr;
2242                                                         VlakRen *vlr= RE_findOrAddVlak(obr, (ps->facenr-1) & RE_QUAD_MASK);
2243                                                         
2244                                                         samp= samplebuf[sample] + sindex;
2245                                                         /* convert image plane pixel location to lamp buffer space */
2246                                                         if(viewpixel_to_lampbuf(shb, obi, vlr, xs + R.jit[sample][0], ys + R.jit[sample][1], samp->zco)) {
2247                                                                 samp->obi= ps->obi;
2248                                                                 samp->facenr= ps->facenr & ~RE_QUAD_OFFS;
2249                                                                 ps->shadfac= 0;
2250                                                                 samp->shadfac= &ps->shadfac;
2251                                                                 bound_rectf((rctf *)&root.box, samp->zco);
2252                                                         }
2253                                                 }
2254                                         }
2255                                 }
2256                         }
2257                         else {
2258                                 rectp= pa->rectp + sindex;
2259                                 recto= pa->recto + sindex;
2260                                 if(*rectp>0) {
2261                                         ObjectInstanceRen *obi= &R.objectinstance[*recto];
2262                                         ObjectRen *obr= obi->obr;
2263                                         VlakRen *vlr= RE_findOrAddVlak(obr, (*rectp-1) & RE_QUAD_MASK);
2264                                         float xs= (float)(x + pa->disprect.xmin);
2265                                         float ys= (float)(y + pa->disprect.ymin);
2266                                         
2267                                         samp= samplebuf[0] + sindex;
2268                                         /* convert image plane pixel location to lamp buffer space */
2269                                         if(viewpixel_to_lampbuf(shb, obi, vlr, xs, ys, samp->zco)) {
2270                                                 samp->obi= *recto;
2271                                                 samp->facenr= *rectp & ~RE_QUAD_OFFS;
2272                                                 samp->shadfac= isbdata->shadfacs + sindex;
2273                                                 bound_rectf((rctf *)&root.box, samp->zco);
2274                                         }
2275                                 }
2276                         }
2277                 }
2278         }
2279         
2280         /* simple method to see if we have samples */
2281         if(root.box.xmin != (float)shb->size) {
2282                 /* now create a regular split, root.box has the initial bounding box of all pixels */
2283                 /* split bsp 8 levels deep, in regular grid (16 x 16) */
2284                 isb_bsp_split_init(&root, memarena, 8);
2285                 
2286                 /* insert all samples in BSP now */
2287                 bsp_err= isb_add_samples(pa, &root, memarena, samplebuf);
2288                         
2289                 if(bsp_err==0) {
2290                         /* go over all faces and fill in shadow values */
2291                         
2292                         isb_bsp_fillfaces(&R, lar, &root);      /* shb->persmat should have been calculated */
2293                         
2294                         /* copy shadow samples to persistant buffer, reduce memory overhead */
2295                         if(R.osa) {
2296                                 ISBShadfacA **isbsa= isbdata->shadfaca= MEM_callocN(pa->rectx*pa->recty*sizeof(void *), "isb shadfacs");
2297                                 
2298                                 isbdata->memarena = BLI_memarena_new(0x8000 * sizeof(ISBSampleA));
2299                                 BLI_memarena_use_calloc(isbdata->memarena);
2300
2301                                 for(rd= pa->rectdaps, x=pa->rectx*pa->recty; x>0; x--, rd++, isbsa++) {
2302                                         
2303                                         if(*rd) {
2304                                                 PixStr *ps= (PixStr *)(*rd);
2305                                                 while(ps) {
2306                                                         if(ps->shadfac)
2307                                                                 isb_add_shadfac(isbsa, isbdata->memarena, ps->obi, ps->facenr, ps->shadfac, count_mask(ps->mask));
2308                                                         ps= ps->next;
2309                                                 }
2310                                         }
2311                                 }
2312                         }
2313                 }
2314         }
2315         else {
2316                 if(isbdata->shadfacs) {
2317                         MEM_freeN(isbdata->shadfacs);
2318                         isbdata->shadfacs= NULL;
2319                 }
2320         }
2321
2322         /* free BSP */
2323         BLI_memarena_free(memarena);
2324         
2325         /* free samples */
2326         for(x=0; x<(R.osa?R.osa:1); x++)
2327                 MEM_freeN(samplebuf[x]);
2328         
2329         if(bsp_err) printf("error in filling bsp\n");
2330 }
2331
2332 /* add sample to buffer, isbsa is the root sample in a buffer */
2333 static ISBSampleA *isb_alloc_sample_transp(ISBSampleA **isbsa, MemArena *mem)
2334 {
2335         ISBSampleA *new;
2336         
2337         new= BLI_memarena_alloc(mem, sizeof(ISBSampleA));
2338         if(*isbsa)
2339                 new->next= (*isbsa);
2340         else
2341                 new->next= NULL;
2342         
2343         *isbsa= new;
2344         return new;
2345 }
2346
2347 /* adding samples in BSP, transparent case */
2348 static int isb_add_samples_transp(RenderPart *pa, ISBBranch *root, MemArena *memarena, ISBSampleA ***samplebuf)
2349 {
2350         int xi, yi, *xcos, *ycos;
2351         int sample, bsp_err= 0;
2352         
2353         /* bsp split doesn't like to handle regular sequenes */
2354         xcos= MEM_mallocN( pa->rectx*sizeof(int), "xcos");
2355         ycos= MEM_mallocN( pa->recty*sizeof(int), "ycos");
2356         for(xi=0; xi<pa->rectx; xi++)
2357                 xcos[xi]= xi;
2358         for(yi=0; yi<pa->recty; yi++)
2359                 ycos[yi]= yi;
2360         BLI_array_randomize(xcos, sizeof(int), pa->rectx, 12345);
2361         BLI_array_randomize(ycos, sizeof(int), pa->recty, 54321);
2362         
2363         for(sample=0; sample<(R.osa?R.osa:1); sample++) {
2364                 ISBSampleA **samp= samplebuf[sample], *samp1;
2365                 
2366                 for(yi=0; yi<pa->recty; yi++) {
2367                         int y= ycos[yi];
2368                         for(xi=0; xi<pa->rectx; xi++) {
2369                                 int x= xcos[xi];
2370                                 
2371                                 samp1= *(samp + y*pa->rectx + x);
2372                                 while(samp1) {
2373                                         bsp_err |= isb_bsp_insert(root, memarena, (ISBSample *)samp1);
2374                                         samp1= samp1->next;
2375                                 }
2376                         }
2377                         if(bsp_err) break;
2378                 }
2379         }       
2380         
2381         MEM_freeN(xcos);
2382         MEM_freeN(ycos);
2383         
2384         return bsp_err;
2385 }
2386
2387
2388 /* Ztransp version */
2389 /* lar->shb, pa->rectz and pa->rectp should exist */
2390 static void isb_make_buffer_transp(RenderPart *pa, APixstr *apixbuf, LampRen *lar)
2391 {
2392         ShadBuf *shb= lar->shb;
2393         ISBData *isbdata;
2394         ISBSampleA *samp, **samplebuf[16];      /* MAX_OSA */
2395         ISBBranch root;
2396         MemArena *memarena;
2397         APixstr *ap;
2398         int x, y, sindex, sample, bsp_err=0;
2399         
2400         /* storage for shadow, per thread */
2401         isbdata= shb->isb_result[pa->thread];
2402         
2403         /* to map the shi->xs and ys coordinate */
2404         isbdata->minx= pa->disprect.xmin;
2405         isbdata->miny= pa->disprect.ymin;
2406         isbdata->rectx= pa->rectx;
2407         isbdata->recty= pa->recty;
2408         
2409         /* branches are added using memarena (32k branches) */
2410         memarena = BLI_memarena_new(0x8000 * sizeof(ISBBranch));
2411         BLI_memarena_use_calloc(memarena);
2412         
2413         /* samplebuf is in camera view space (pixels) */
2414         for(sample=0; sample<(R.osa?R.osa:1); sample++)
2415                 samplebuf[sample]= MEM_callocN(sizeof(void *)*pa->rectx*pa->recty, "isb alpha samplebuf");
2416         
2417         /* setup bsp root */
2418         memset(&root, 0, sizeof(ISBBranch));
2419         root.box.xmin= (float)shb->size;
2420         root.box.ymin= (float)shb->size;
2421
2422         /* create the sample buffers */
2423         for(ap= apixbuf, sindex=0, y=0; y<pa->recty; y++) {
2424                 for(x=0; x<pa->rectx; x++, sindex++, ap++) {
2425                         
2426                         if(ap->p[0]) {
2427                                 APixstr *apn;
2428                                 float xs= (float)(x + pa->disprect.xmin);
2429                                 float ys= (float)(y + pa->disprect.ymin);
2430                                 
2431                                 for(apn=ap; apn; apn= apn->next) {
2432                                         int a;
2433                                         for(a=0; a<4; a++) {
2434                                                 if(apn->p[a]) {
2435                                                         ObjectInstanceRen *obi= &R.objectinstance[apn->obi[a]];
2436                                                         ObjectRen *obr= obi->obr;
2437                                                         VlakRen *vlr= RE_findOrAddVlak(obr, (apn->p[a]-1) & RE_QUAD_MASK);
2438                                                         float zco[3];
2439                                                         
2440                                                         /* here we store shadfac, easier to create the end storage buffer. needs zero'ed, multiple shadowbufs use it */
2441                                                         apn->shadfac[a]= 0;
2442                                                         
2443                                                         if(R.osa) {
2444                                                                 for(sample=0; sample<R.osa; sample++) {
2445                                                                         int mask= (1<<sample);
2446                                                                         
2447                                                                         if(apn->mask[a] & mask) {
2448                                                                                 
2449                                                                                 /* convert image plane pixel location to lamp buffer space */
2450                                                                                 if(viewpixel_to_lampbuf(shb, obi, vlr, xs + R.jit[sample][0], ys + R.jit[sample][1], zco)) {
2451                                                                                         samp= isb_alloc_sample_transp(samplebuf[sample] + sindex, memarena);
2452                                                                                         samp->obi= apn->obi[a];
2453                                                                                         samp->facenr= apn->p[a] & ~RE_QUAD_OFFS;
2454                                                                                         samp->shadfac= &apn->shadfac[a];
2455                                                                                         
2456                                                                                         VECCOPY(samp->zco, zco);
2457                                                                                         bound_rectf((rctf *)&root.box, samp->zco);
2458                                                                                 }
2459                                                                         }
2460                                                                 }
2461                                                         }
2462                                                         else {
2463                                                                 
2464                                                                 /* convert image plane pixel location to lamp buffer space */
2465                                                                 if(viewpixel_to_lampbuf(shb, obi, vlr, xs, ys, zco)) {
2466                                                                         
2467                                                                         samp= isb_alloc_sample_transp(samplebuf[0] + sindex, memarena);
2468                                                                         samp->obi= apn->obi[a];
2469                                                                         samp->facenr= apn->p[a] & ~RE_QUAD_OFFS;
2470                                                                         samp->shadfac= &apn->shadfac[a];
2471                                                                         
2472                                                                         VECCOPY(samp->zco, zco);
2473                                                                         bound_rectf((rctf *)&root.box, samp->zco);
2474                                                                 }
2475                                                         }
2476                                                 }
2477                                         }
2478                                 }
2479                         }
2480                 }
2481         }
2482         
2483         /* simple method to see if we have samples */
2484         if(root.box.xmin != (float)shb->size) {
2485                 /* now create a regular split, root.box has the initial bounding box of all pixels */
2486                 /* split bsp 8 levels deep, in regular grid (16 x 16) */
2487                 isb_bsp_split_init(&root, memarena, 8);
2488                 
2489                 /* insert all samples in BSP now */
2490                 bsp_err= isb_add_samples_transp(pa, &root, memarena, samplebuf);
2491                 
2492                 if(bsp_err==0) {
2493                         ISBShadfacA **isbsa;
2494                         
2495                         /* go over all faces and fill in shadow values */
2496                         isb_bsp_fillfaces(&R, lar, &root);      /* shb->persmat should have been calculated */
2497                         
2498                         /* copy shadow samples to persistant buffer, reduce memory overhead */
2499                         isbsa= isbdata->shadfaca= MEM_callocN(pa->rectx*pa->recty*sizeof(void *), "isb shadfacs");
2500                         
2501                         isbdata->memarena = BLI_memarena_new(0x8000 * sizeof(ISBSampleA));
2502                         
2503                         for(ap= apixbuf, x=pa->rectx*pa->recty; x>0; x--, ap++, isbsa++) {
2504                                         
2505                                 if(ap->p[0]) {
2506                                         APixstr *apn;
2507                                         for(apn=ap; apn; apn= apn->next) {
2508                                                 int a;
2509                                                 for(a=0; a<4; a++) {
2510                                                         if(apn->p[a] && apn->shadfac[a]) {
2511                                                                 if(R.osa)
2512                                                                         isb_add_shadfac(isbsa, isbdata->memarena, apn->obi[a], apn->p[a], apn->shadfac[a], count_mask(apn->mask[a]));
2513                                                                 else
2514                                                                         isb_add_shadfac(isbsa, isbdata->memarena, apn->obi[a], apn->p[a], apn->shadfac[a], 0);
2515                                                         }
2516                                                 }
2517                                         }
2518                                 }
2519                         }
2520                 }
2521         }
2522
2523         /* free BSP */
2524         BLI_memarena_free(memarena);
2525
2526         /* free samples */
2527         for(x=0; x<(R.osa?R.osa:1); x++)
2528                 MEM_freeN(samplebuf[x]);
2529
2530         if(bsp_err) printf("error in filling bsp\n");
2531 }
2532
2533
2534
2535 /* exported */
2536
2537 /* returns amount of light (1.0 = no shadow) */
2538 /* note, shadepixel() rounds the coordinate, not the real sample info */
2539 float ISB_getshadow(ShadeInput *shi, ShadBuf *shb)
2540 {
2541         /* if raytracing, we can't accept irregular shadow */
2542         if(shi->depth==0) {
2543                 ISBData *isbdata= shb->isb_result[shi->thread];
2544                 
2545                 if(isbdata) {
2546                         if(isbdata->shadfacs || isbdata->shadfaca) {
2547                                 int x= shi->xs - isbdata->minx;
2548                                 
2549                                 if(x >= 0 && x < isbdata->rectx) {
2550                                         int y= shi->ys - isbdata->miny;
2551                         
2552                                         if(y >= 0 && y < isbdata->recty) {
2553                                                 if(isbdata->shadfacs) {
2554                                                         short *sp= isbdata->shadfacs + y*isbdata->rectx + x;
2555                                                         return *sp>=4096?0.0f:1.0f - ((float)*sp)/4096.0f;
2556                                                 }
2557                                                 else {
2558                                                         int sindex= y*isbdata->rectx + x;
2559                                                         int obi= shi->obi - R.objectinstance;
2560                                                         ISBShadfacA *isbsa= *(isbdata->shadfaca + sindex);
2561                                                         
2562                                                         while(isbsa) {
2563                                                                 if(isbsa->facenr==shi->facenr+1 && isbsa->obi==obi)
2564                                                                         return isbsa->shadfac>=1.0f?0.0f:1.0f - isbsa->shadfac;
2565                                                                 isbsa= isbsa->next;
2566                                                         }
2567                                                 }
2568                                         }
2569                                 }
2570                         }
2571                 }
2572         }
2573         return 1.0f;
2574 }
2575
2576 /* part is supposed to be solid zbuffered (apixbuf==NULL) or transparent zbuffered */
2577 void ISB_create(RenderPart *pa, APixstr *apixbuf)
2578 {
2579         GroupObject *go;
2580         
2581         /* go over all lamps, and make the irregular buffers */
2582         for(go=R.lights.first; go; go= go->next) {
2583                 LampRen *lar= go->lampren;
2584                 
2585                 if(lar->type==LA_SPOT && lar->shb && lar->buftype==LA_SHADBUF_IRREGULAR) {
2586                         
2587                         /* create storage for shadow, per thread */
2588                         lar->shb->isb_result[pa->thread]= MEM_callocN(sizeof(ISBData), "isb data");
2589                         
2590                         if(apixbuf)
2591                                 isb_make_buffer_transp(pa, apixbuf, lar);
2592                         else
2593                                 isb_make_buffer(pa, lar);
2594                 }
2595         }
2596 }
2597
2598
2599 /* end of part rendering, free stored shadow data for this thread from all lamps */
2600 void ISB_free(RenderPart *pa)
2601 {
2602         GroupObject *go;
2603         
2604         /* go over all lamps, and free the irregular buffers */
2605         for(go=R.lights.first; go; go= go->next) {
2606                 LampRen *lar= go->lampren;
2607                 
2608                 if(lar->type==LA_SPOT && lar->shb && lar->buftype==LA_SHADBUF_IRREGULAR) {
2609                         ISBData *isbdata= lar->shb->isb_result[pa->thread];
2610
2611                         if(isbdata) {
2612                                 if(isbdata->shadfacs)
2613                                         MEM_freeN(isbdata->shadfacs);
2614                                 if(isbdata->shadfaca)
2615                                         MEM_freeN(isbdata->shadfaca);
2616                                 
2617                                 if(isbdata->memarena)
2618                                         BLI_memarena_free(isbdata->memarena);
2619                                 
2620                                 MEM_freeN(isbdata);
2621                                 lar->shb->isb_result[pa->thread]= NULL;
2622                         }
2623                 }
2624         }
2625 }
2626