The zblur plugin (aka as DoF) integrated in render. Compared to patch
[blender.git] / source / blender / render / intern / source / shadbuf.c
1 /*  shadbuf.c        RENDER
2  *
3  * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version. The Blender
9  * Foundation also sells licenses for use in proprietary software under
10  * the Blender License.  See http://www.blender.org/BL/ for information
11  * about this.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software Foundation,
20  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
21  *
22  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
23  * All rights reserved.
24  *
25  * The Original Code is: all of this file.
26  *
27  * Contributor(s): none yet.
28  *
29  * ***** END GPL/BL DUAL LICENSE BLOCK *****
30  * 
31  *  april 95
32  * 
33  * $Id$
34  *
35  * 27-Jun-2001 switched the shadow buffer for UR to the object-shadow
36  * buffers, and removed all references and fixes for UR rendering from
37  * this one.
38  * */
39
40 #include <math.h>
41 #include <string.h>
42
43 #include "MTC_matrixops.h"
44 #include "MEM_guardedalloc.h"
45
46 #include "DNA_lamp_types.h"
47 #include "BKE_utildefines.h"
48 #include "BLI_arithb.h"
49
50 #include "render.h"
51
52 #include "shadbuf.h"
53 #include "renderHelp.h"
54 #include "jitter.h"
55 #include "zbuf.h"
56
57 /* XXX, could be better implemented...
58 */
59 #if defined(__sgi) || defined(__sparc) || defined(__sparc__) || defined (__PPC__) || defined (__ppc__) || defined (__BIG_ENDIAN__)
60 #define RCOMP   3
61 #define GCOMP   2
62 #define BCOMP   1
63 #define ACOMP   0
64 #else
65 #define RCOMP   0
66 #define GCOMP   1
67 #define BCOMP   2
68 #define ACOMP   3
69 #endif
70
71 /* ------------------------------------------------------------------------- */
72
73
74 void RE_initshadowbuf(LampRen *lar, float mat[][4])
75 {
76         struct ShadBuf *shb;
77         float hoek, temp, viewinv[4][4];
78
79         /* if(la->spsi<16) return; */
80
81         /* memory reservation */
82         shb= (struct ShadBuf *)MEM_callocN( sizeof(struct ShadBuf),"initshadbuf");
83         lar->shb= shb;
84
85         if(shb==NULL) return;
86
87         VECCOPY(shb->co, lar->co);
88         
89         /* percentage render: keep track of min and max */
90         shb->size= (lar->bufsize*R.r.size)/100;
91         if(shb->size<512) shb->size= 512;
92         else if(shb->size > lar->bufsize) shb->size= lar->bufsize;
93
94         shb->size &= ~15;       /* make sure its multiples of 16 */
95
96         shb->samp= lar->samp;
97         shb->soft= lar->soft;
98         shb->shadhalostep= lar->shadhalostep;
99         
100         shb->zbuf= (unsigned long *)MEM_mallocN( sizeof(unsigned long)*(shb->size*shb->size)/256, "initshadbuf2");
101         shb->cbuf= (char *)MEM_callocN( (shb->size*shb->size)/256, "initshadbuf3");
102
103         if(shb->zbuf==0 || shb->cbuf==0) {
104                 if(shb->zbuf) MEM_freeN(shb->zbuf);
105                 MEM_freeN(lar->shb);
106                 lar->shb= 0;            
107                 return;
108         }
109
110         MTC_Mat4Ortho(mat);
111         MTC_Mat4Invert(shb->winmat, mat);       /* winmat is temp */
112
113         /* matrix: combination of inverse view and lampmat */
114         /* calculate again: the ortho-render has no correct viewinv */
115         MTC_Mat4Invert(viewinv, R.viewmat);
116         MTC_Mat4MulMat4(shb->viewmat, viewinv, shb->winmat);
117
118         /* projection */
119         hoek= saacos(lar->spotsi);
120         temp= 0.5*shb->size*cos(hoek)/sin(hoek);
121         shb->d= lar->clipsta;
122
123         shb->pixsize= (shb->d)/temp;
124
125         shb->far= lar->clipend;
126         /* bias is percentage, made 2x karger because of correction for angle of incidence */
127         /* when a ray is closer to parallel of a face, bias value is increased during render */
128         shb->bias= (0.02*lar->bias)*0x7FFFFFFF;
129         shb->bias= shb->bias*(100/R.r.size);
130         
131 }
132 /* ------------------------------------------------------------------------- */
133
134
135 static void lrectreadRectz(int x1, int y1, int x2, int y2, char *r1) /* reads part from rectz in r1 */
136 {
137         unsigned int len4, *rz; 
138
139         if(x1>=R.rectx || x2>=R.rectx || y1>=R.recty || y2>=R.recty) return;
140         if(x1>x2 || y1>y2) return;
141
142         len4= 4*(x2- x1+1);
143         rz= R.rectz+R.rectx*y1+x1;
144         for(;y1<=y2;y1++) {
145                 memcpy(r1,rz,len4);
146                 rz+= R.rectx;
147                 r1+= len4;
148         }
149 }
150
151 #if 0
152 static int sizeoflampbuf(struct ShadBuf *shb)
153 {
154         int num,count=0;
155         char *cp;
156         
157         cp= shb->cbuf;
158         num= (shb->size*shb->size)/256;
159
160         while(num--) count+= *(cp++);
161         
162         return 256*count;
163 }
164 #endif
165
166 static float *give_jitter_tab(int samp)
167 {
168         /* these are all possible jitter tables, takes up some
169          * 12k, not really bad!
170          * For soft shadows, it saves memory and render time
171          */
172         static int tab[17]={1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, 196, 225, 256};
173         static float jit[1496][2];
174         static char ctab[17]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
175         int a, offset=0;
176         
177         if(samp<2) samp= 2;
178         else if(samp>16) samp= 16;
179
180         for(a=0; a<samp-1; a++) offset+= tab[a];
181
182         if(ctab[samp]==0) {
183                 initjit(jit[offset], samp*samp);
184                 ctab[samp]= 1;
185         }
186                 
187         return jit[offset];
188         
189 }
190
191 void makeshadowbuf(LampRen *lar)
192 {
193         struct ShadBuf *shb= lar->shb;
194         float panophi;
195         float temp, wsize, dist;
196         int *rz, *rz1, verg, verg1;
197         unsigned long *ztile;
198         int a, x, y, minx, miny, byt1, byt2, tempmode;
199         short temprx,tempry, square;
200         char *rc, *rcline, *ctile, *zt;
201
202         panophi = getPanoPhi();
203         
204         /* store viewvars */
205         temprx= R.rectx; tempry= R.recty; 
206         tempmode= R.r.mode;
207         R.r.mode &= ~R_ORTHO;
208         R.rectx= R.recty= shb->size;
209
210         shb->jit= give_jitter_tab(shb->samp);
211
212         /* matrices and window: in R.winmat the transformation is being put,
213                 transforming from observer view to lamp view, including lamp window matrix */
214         
215         wsize= shb->pixsize*(shb->size/2.0);
216
217         i_window(-wsize, wsize, -wsize, wsize, shb->d, shb->far, shb->winmat);
218
219         MTC_Mat4MulMat4(shb->persmat, shb->viewmat, shb->winmat);
220         
221         /* temp, will be restored */
222         MTC_Mat4SwapMat4(shb->persmat, R.winmat);
223
224         /* zbuffering */
225         if(R.rectz) MEM_freeN(R.rectz);
226         R.rectz= (unsigned int *)MEM_mallocN(sizeof(int)*shb->size*shb->size,"makeshadbuf");
227         rcline= MEM_mallocN(256*4+sizeof(int),"makeshadbuf2");
228
229         /* store: panorama rot */
230         temp= panophi;
231         panophi= 0.0;
232         pushTempPanoPhi(0.0);
233
234         /* pano interference here? */
235         setzbufvlaggen(projectvert);
236
237         popTempPanoPhi();
238         panophi= temp;
239         
240         zbuffershad(lar);
241         
242         square= lar->mode & LA_SQUARE;
243
244         /* create Z tiles (for compression): this system is 24 bits!!! */
245         
246         ztile= shb->zbuf;
247         ctile= shb->cbuf;
248         for(y=0; y<shb->size; y+=16) {
249                 if(y< shb->size/2) miny= y+15-shb->size/2;
250                 else miny= y-shb->size/2;       
251                                 
252                 for(x=0; x<shb->size; x+=16) {
253
254                         /* is tile within spotbundle? */
255                         a= shb->size/2;
256                         if(x< a) minx= x+15-a;
257                         else minx= x-a; 
258                         
259                         dist= sqrt( (float)(minx*minx+miny*miny) );
260
261                         if(square==0 && dist>(float)(a+12)) {   /* 12, tested with a onlyshadow lamp */
262                                 a= 256; verg= 0; /* 0x80000000; */ /* 0x7FFFFFFF; */
263                                 rz1= (&verg)+1;
264                         } 
265                         else {
266                                 lrectreadRectz(x, y, MIN2(shb->size-1,x+15), MIN2(shb->size-1,y+15), rcline);
267                                 rz1= (int *)rcline;
268                                 
269                                 verg= (*rz1 & 0xFFFFFF00);
270
271                                 for(a=0;a<256;a++,rz1++) {
272                                         if( (*rz1 & 0xFFFFFF00) !=verg) break;
273                                 }
274                         }
275                         if(a==256) { /* complete empty tile */
276                                 *ctile= 0;
277                                 *ztile= *(rz1-1);
278                         }
279                         else {
280
281                                 /* ACOMP etc. are defined to work L/B endian */
282                                 
283                                 rc= rcline;
284                                 rz1= (int *)rcline;
285                                 verg=  rc[ACOMP];
286                                 verg1= rc[BCOMP];
287                                 rc+= 4;
288                                 byt1= 1; byt2= 1;
289                                 for(a=1;a<256;a++,rc+=4) {
290                                         byt1 &= (verg==rc[ACOMP]);
291                                         byt2 &= (verg1==rc[BCOMP]);
292
293                                         if(byt1==0) break;
294                                 }
295                                 if(byt1 && byt2) {      /* only store byte */
296                                         *ctile= 1;
297                                         *ztile= (unsigned long)MEM_mallocN(256+4, "tile1");
298                                         rz= (int *)*ztile;
299                                         *rz= *rz1;
300
301                                         zt= (char *)(rz+1);
302                                         rc= rcline;
303                                         for(a=0; a<256; a++, zt++, rc+=4) *zt= rc[GCOMP];       
304                                 }
305                                 else if(byt1) {         /* only store short */
306                                         *ctile= 2;
307                                         *ztile= (unsigned long)MEM_mallocN(2*256+4,"Tile2");
308                                         rz= (int *)*ztile;
309                                         *rz= *rz1;
310
311                                         zt= (char *)(rz+1);
312                                         rc= rcline;
313                                         for(a=0; a<256; a++, zt+=2, rc+=4) {
314                                                 zt[0]= rc[BCOMP];
315                                                 zt[1]= rc[GCOMP];
316                                         }
317                                 }
318                                 else {                  /* store triple */
319                                         *ctile= 3;
320                                         *ztile= (unsigned long)MEM_mallocN(3*256,"Tile3");
321
322                                         zt= (char *)*ztile;
323                                         rc= rcline;
324                                         for(a=0; a<256; a++, zt+=3, rc+=4) {
325                                                 zt[0]= rc[ACOMP];
326                                                 zt[1]= rc[BCOMP];
327                                                 zt[2]= rc[GCOMP];
328                                         }
329                                 }
330                         }
331                         ztile++;
332                         ctile++;
333                 }
334         }
335
336         MEM_freeN(rcline);
337         MEM_freeN(R.rectz); R.rectz= NULL;
338         
339         /* old globals back */
340         R.rectx= temprx; R.recty= tempry;
341         R.r.mode= tempmode;
342         MTC_Mat4SwapMat4(shb->persmat, R.winmat);
343
344         /* printf("lampbuf %d\n", sizeoflampbuf(shb)); */
345 }
346
347 static int firstreadshadbuf(struct ShadBuf *shb, int xs, int ys, int nr)
348 {
349         /* return a 1 if fully compressed shadbuf-tile && z==const */
350         static int *rz;
351         int ofs;
352         char *ct;
353
354         /* always test borders of shadowbuffer */
355         if(xs<0) xs= 0; else if(xs>=shb->size) xs= shb->size-1;
356         if(ys<0) ys= 0; else if(ys>=shb->size) ys= shb->size-1;
357    
358         /* calc z */
359         ofs= (ys>>4)*(shb->size>>4) + (xs>>4);
360         ct= shb->cbuf+ofs;
361         if(*ct==0) {
362             if(nr==0) {
363                         rz= *( (int **)(shb->zbuf+ofs) );
364                         return 1;
365             }
366                 else if(rz!= *( (int **)(shb->zbuf+ofs) )) return 0;
367                 
368             return 1;
369         }
370         
371         return 0;
372 }
373
374 static float readshadowbuf(struct ShadBuf *shb, int bias, int xs, int ys, int zs)       /* return 1.0 : fully in light */
375 {
376         float temp;
377         int *rz, ofs;
378         int zsamp=0;
379         char *ct, *cz;
380
381         /* simpleclip */
382         /* if(xs<0 || ys<0) return 1.0; */
383         /* if(xs>=shb->size || ys>=shb->size) return 1.0; */
384         
385         /* always test borders of shadowbuffer */
386         if(xs<0) xs= 0; else if(xs>=shb->size) xs= shb->size-1;
387         if(ys<0) ys= 0; else if(ys>=shb->size) ys= shb->size-1;
388
389         /* calc z */
390         ofs= (ys>>4)*(shb->size>>4) + (xs>>4);
391         ct= shb->cbuf+ofs;
392         rz= *( (int **)(shb->zbuf+ofs) );
393
394         if(*ct==3) {
395                 ct= ((char *)rz)+3*16*(ys & 15)+3*(xs & 15);
396                 cz= (char *)&zsamp;
397                 cz[ACOMP]= ct[0];
398                 cz[BCOMP]= ct[1];
399                 cz[GCOMP]= ct[2];
400         }
401         else if(*ct==2) {
402                 ct= ((char *)rz);
403                 ct+= 4+2*16*(ys & 15)+2*(xs & 15);
404                 zsamp= *rz;
405         
406                 cz= (char *)&zsamp;
407                 cz[BCOMP]= ct[0];
408                 cz[GCOMP]= ct[1];
409         }
410         else if(*ct==1) {
411                 ct= ((char *)rz);
412                 ct+= 4+16*(ys & 15)+(xs & 15);
413                 zsamp= *rz;
414
415                 cz= (char *)&zsamp;
416                 cz[GCOMP]= ct[0];
417
418         }
419         else {
420                 /* got warning on this from DEC alpha (64 bits).... */
421                 /* but it's working code! (ton) */
422                 zsamp= (int) rz;
423         }
424
425         /* if(zsamp >= 0x7FFFFE00) return 1.0; */       /* no shaduw when sampling at eternal distance */
426
427         if(zsamp > zs) return 1.0;              /* absolute no shadow */
428         else if( zsamp < zs-bias) return 0.0 ;  /* absolute in shadow */
429         else {                                  /* soft area */
430                 
431                 temp=  ( (float)(zs- zsamp) )/(float)bias;
432                 return 1.0 - temp*temp;
433                         
434         }
435 }
436
437 /* the externally called shadow testing (reading) function */
438 /* return 1.0: no shadow at all */
439 float testshadowbuf(struct ShadBuf *shb, float *rco, float *dxco, float *dyco, float inp)
440 {
441         float fac, co[4], dx[3], dy[3], aantal=0;
442         float xs1,ys1, siz, *j, xres, yres;
443         int xs,ys, zs, bias;
444         short a,num;
445         
446         /* if(inp <= 0.0) return 1.0; */
447
448         /* rotate renderco en osaco */
449         siz= 0.5*(float)shb->size;
450         VECCOPY(co, rco);
451         co[3]= 1.0;
452
453         MTC_Mat4MulVec4fl(shb->persmat, co);    /* rational hom co */
454
455         xs1= siz*(1.0+co[0]/co[3]);
456         ys1= siz*(1.0+co[1]/co[3]);
457
458         /* Clip for z: near and far clip values of the shadow buffer. We
459                 * can test for -1.0/1.0 because of the properties of the
460                 * coordinate transformations. */
461         fac= (co[2]/co[3]);
462
463         if(fac>=1.0) {
464                 return 0.0;
465         } else if(fac<= -1.0) {
466                 return 1.0;
467         }
468
469         zs= ((float)0x7FFFFFFF)*fac;
470
471         /* take num*num samples, increase area with fac */
472         num= shb->samp*shb->samp;
473         fac= shb->soft;
474         
475         /* with inp==1.0, bias is half the size. correction value was 1.1, giving errors 
476            on cube edges, with one side being almost frontal lighted (ton)  */
477         bias= (1.5-inp*inp)*shb->bias;
478
479         if(num==1) {
480                 return readshadowbuf(shb, bias, (int)xs1, (int)ys1, zs);
481         }
482
483         co[0]= rco[0]+dxco[0];
484         co[1]= rco[1]+dxco[1];
485         co[2]= rco[2]+dxco[2];
486         co[3]= 1.0;
487         MTC_Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
488         dx[0]= xs1- siz*(1.0+co[0]/co[3]);
489         dx[1]= ys1- siz*(1.0+co[1]/co[3]);
490
491         co[0]= rco[0]+dyco[0];
492         co[1]= rco[1]+dyco[1];
493         co[2]= rco[2]+dyco[2];
494         co[3]= 1.0;
495         MTC_Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
496         dy[0]= xs1- siz*(1.0+co[0]/co[3]);
497         dy[1]= ys1- siz*(1.0+co[1]/co[3]);
498
499         xres= fac*( fabs(dx[0])+fabs(dy[0]) );
500         yres= fac*( fabs(dx[1])+fabs(dy[1]) );
501
502         if(xres<fac) xres= fac;
503         if(yres<fac) yres= fac;
504         
505         xs1-= (xres)/2;
506         ys1-= (yres)/2;
507
508         j= shb->jit;
509
510         if(xres<16.0 && yres<16.0) {
511             if(firstreadshadbuf(shb, (int)xs1, (int)ys1, 0)) {
512                         if(firstreadshadbuf(shb, (int)(xs1+xres), (int)ys1, 1)) {
513                                 if(firstreadshadbuf(shb, (int)xs1, (int)(ys1+yres), 1)) {
514                                         if(firstreadshadbuf(shb, (int)(xs1+xres), (int)(ys1+yres), 1)) {
515                                                 return readshadowbuf(shb, bias,(int)xs1, (int)ys1, zs);
516                                         }
517                                 }
518                         }
519             }
520         }
521
522         for(a=num;a>0;a--) {
523                 /* instead of jit i tried random: ugly! */
524                 /* note: the plus 0.5 gives best sampling results, jit used to go from 0-1 */
525                 /* xs1 and ys1 are already corrected to be corner of sample area */
526                 xs= xs1 + xres*(j[0] + 0.5);
527                 ys= ys1 + yres*(j[1] + 0.5);
528                 j+=2;
529                 
530                 aantal+= readshadowbuf(shb, bias, xs, ys, zs);
531         }
532
533         /* Renormalizes for the sample number: */
534         return aantal/( (float)(num) );
535 }
536
537 /* different function... sampling behind clipend can be LIGHT, bias is negative! */
538 /* return: light */
539 static float readshadowbuf_halo(struct ShadBuf *shb, int xs, int ys, int zs)
540 {
541         float temp;
542         int *rz, ofs;
543         int bias, zbias, zsamp;
544         char *ct, *cz;
545
546         /* negative! The other side is more important */
547         bias= -shb->bias;
548         
549         /* simpleclip */
550         if(xs<0 || ys<0) return 0.0;
551         if(xs>=shb->size || ys>=shb->size) return 0.0;
552
553         /* calc z */
554         ofs= (ys>>4)*(shb->size>>4) + (xs>>4);
555         ct= shb->cbuf+ofs;
556         rz= *( (int **)(shb->zbuf+ofs) );
557
558         if(*ct==3) {
559                 ct= ((char *)rz)+3*16*(ys & 15)+3*(xs & 15);
560                 cz= (char *)&zsamp;
561                 zsamp= 0;
562                 cz[ACOMP]= ct[0];
563                 cz[BCOMP]= ct[1];
564                 cz[GCOMP]= ct[2];
565         }
566         else if(*ct==2) {
567                 ct= ((char *)rz);
568                 ct+= 4+2*16*(ys & 15)+2*(xs & 15);
569                 zsamp= *rz;
570         
571                 cz= (char *)&zsamp;
572                 cz[BCOMP]= ct[0];
573                 cz[GCOMP]= ct[1];
574         }
575         else if(*ct==1) {
576                 ct= ((char *)rz);
577                 ct+= 4+16*(ys & 15)+(xs & 15);
578                 zsamp= *rz;
579
580                 cz= (char *)&zsamp;
581                 cz[GCOMP]= ct[0];
582
583         }
584         else {
585                 /* same as before */
586                 /* still working code! (ton) */
587                 zsamp= (int) rz;
588         }
589
590         /* NO schadow when sampled at 'eternal' distance */
591
592         if(zsamp >= 0x7FFFFE00) return 1.0; 
593
594         if(zsamp > zs) return 1.0;              /* absolute no shadww */
595         else {
596                 /* bias is negative, so the (zs-bias) can be beyond 0x7fffffff */
597                 zbias= 0x7fffffff - zs;
598                 if(zbias > -bias) {
599                         if( zsamp < zs-bias) return 0.0 ;       /* absolute in shadow */
600                 }
601                 else return 0.0 ;       /* absolute shadow */
602         }
603
604         /* soft area */
605         
606         temp=  ( (float)(zs- zsamp) )/(float)bias;
607         return 1.0 - temp*temp;
608 }
609
610
611 float shadow_halo(LampRen *lar, float *p1, float *p2)
612 {
613         /* p1 p2 already are rotated in spot-space */
614         ShadBuf *shb= lar->shb;
615         float co[4], siz;
616         float labda, labdao, labdax, labday, ldx, ldy;
617         float zf, xf1, yf1, zf1, xf2, yf2, zf2;
618         float count, lightcount;
619         int x, y, z, xs1, ys1;
620         int dx = 0, dy = 0;
621         
622         siz= 0.5*(float)shb->size;
623         
624         co[0]= p1[0];
625         co[1]= p1[1];
626         co[2]= p1[2]/lar->sh_zfac;
627         co[3]= 1.0;
628         MTC_Mat4MulVec4fl(shb->winmat, co);     /* rational hom co */
629         xf1= siz*(1.0+co[0]/co[3]);
630         yf1= siz*(1.0+co[1]/co[3]);
631         zf1= (co[2]/co[3]);
632
633
634         co[0]= p2[0];
635         co[1]= p2[1];
636         co[2]= p2[2]/lar->sh_zfac;
637         co[3]= 1.0;
638         MTC_Mat4MulVec4fl(shb->winmat, co);     /* rational hom co */
639         xf2= siz*(1.0+co[0]/co[3]);
640         yf2= siz*(1.0+co[1]/co[3]);
641         zf2= (co[2]/co[3]);
642
643         /* the 2dda (a pixel line formula) */
644
645         xs1= (int)xf1;
646         ys1= (int)yf1;
647
648         if(xf1 != xf2) {
649                 if(xf2-xf1 > 0.0) {
650                         labdax= (xf1-xs1-1.0)/(xf1-xf2);
651                         ldx= -shb->shadhalostep/(xf1-xf2);
652                         dx= shb->shadhalostep;
653                 }
654                 else {
655                         labdax= (xf1-xs1)/(xf1-xf2);
656                         ldx= shb->shadhalostep/(xf1-xf2);
657                         dx= -shb->shadhalostep;
658                 }
659         }
660         else {
661                 labdax= 1.0;
662                 ldx= 0.0;
663         }
664
665         if(yf1 != yf2) {
666                 if(yf2-yf1 > 0.0) {
667                         labday= (yf1-ys1-1.0)/(yf1-yf2);
668                         ldy= -shb->shadhalostep/(yf1-yf2);
669                         dy= shb->shadhalostep;
670                 }
671                 else {
672                         labday= (yf1-ys1)/(yf1-yf2);
673                         ldy= shb->shadhalostep/(yf1-yf2);
674                         dy= -shb->shadhalostep;
675                 }
676         }
677         else {
678                 labday= 1.0;
679                 ldy= 0.0;
680         }
681         
682         x= xs1;
683         y= ys1;
684         labda= count= lightcount= 0.0;
685
686 /* printf("start %x %x  \n", (int)(0x7FFFFFFF*zf1), (int)(0x7FFFFFFF*zf2)); */
687
688         while(1) {
689                 labdao= labda;
690                 
691                 if(labdax==labday) {
692                         labdax+= ldx;
693                         x+= dx;
694                         labday+= ldy;
695                         y+= dy;
696                 }
697                 else {
698                         if(labdax<labday) {
699                                 labdax+= ldx;
700                                 x+= dx;
701                         } else {
702                                 labday+= ldy;
703                                 y+= dy;
704                         }
705                 }
706                 
707                 labda= MIN2(labdax, labday);
708                 if(labda==labdao || labda>=1.0) break;
709                 
710                 zf= zf1 + labda*(zf2-zf1);
711                 count+= 1.0;
712
713                 if(zf<= -1.0) lightcount += 1.0;        /* close to the spot */
714                 else {
715                 
716                         /* make sure, behind the clipend we extend halolines. */
717                         if(zf>=1.0) z= 0x7FFFF000;
718                         else z= (int)(0x7FFFF000*zf);
719                         
720                         lightcount+= readshadowbuf_halo(shb, x, y, z);
721                         
722                 }
723         }
724         
725         if(count!=0.0) return (lightcount/count);
726         return 0.0;
727         
728 }
729
730
731
732
733
734
735
736
737