ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond4.c
(Generate patch)

Comparing ray/src/px/pcond4.c (file contents):
Revision 3.2 by greg, Thu Oct 3 20:14:20 1996 UTC vs.
Revision 3.14 by greg, Tue Jan 28 16:31:17 1997 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 1996 Regents of the University of California */
1 > /* Copyright (c) 1997 Regents of the University of California */
2  
3   #ifndef lint
4   static char SCCSid[] = "$SunId$ LBL";
# Line 10 | Line 10 | static char SCCSid[] = "$SunId$ LBL";
10  
11   #include "pcond.h"
12  
13 + /************** VEILING STUFF *****************/
14  
15   #define VADAPT          0.08            /* fraction of adaptation from veil */
16  
17 < extern COLOR    *fovimg;                /* foveal (1 degree) averaged image */
17 < extern short    fvxr, fvyr;             /* foveal image resolution */
17 > static COLOR    *veilimg = NULL;        /* veiling image */
18  
19 – #define fovscan(y)      (fovimg+(y)*fvxr)
20 –
21 – static COLOR    *veilimg;               /* veiling image */
22 –
19   #define veilscan(y)     (veilimg+(y)*fvxr)
20  
21   static float    (*raydir)[3] = NULL;    /* ray direction for each pixel */
# Line 82 | Line 78 | compveil()                             /* compute veiling image */
78          COLOR   ctmp, vsum;
79          int     px, py;
80          register int    x, y;
81 +
82 +        if (veilimg != NULL)            /* already done? */
83 +                return;
84                                          /* compute ray directions */
85          compraydir();
86                                          /* compute veil image */
# Line 98 | Line 97 | compveil()                             /* compute veiling image */
97                                          t2 = DOT(rdirscan(py)[px],
98                                                          rdirscan(y)[x]);
99                                          if (t2 <= FTINY) continue;
100 <                                        t2 = acos(t2);
101 <                                        t2 = 1./(t2*t2);
100 >                                        /*      use approximation instead
101 >                                        t3 = acos(t2);
102 >                                        t2 = t2/(t3*t3);
103 >                                        */
104 >                                        t2 *= .5 / (1. - t2);
105                                          copycolor(ctmp, fovscan(y)[x]);
106                                          scalecolor(ctmp, t2);
107                                          addcolor(vsum, ctmp);
# Line 109 | Line 111 | compveil()                             /* compute veiling image */
111                          scalecolor(vsum, VADAPT/t2sum);
112                          copycolor(veilscan(py)[px], vsum);
113                  }
114 +                                        /* modify FOV sample image */
115 +        for (y = 0; y < fvyr; y++)
116 +                for (x = 0; x < fvxr; x++) {
117 +                        scalecolor(fovscan(y)[x], 1.-VADAPT);
118 +                        addcolor(fovscan(y)[x], veilscan(y)[x]);
119 +                }
120 +        comphist();                     /* recompute histogram */
121   }
122  
123  
# Line 122 | Line 131 | int    y;
131          register int    x, i;
132  
133          vy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
134 <        if (vy >= fvyr-1) vy--;
134 >        while (vy >= fvyr-1) vy--;
135          dy -= (double)vy;
136          for (x = 0; x < scanlen(&inpres); x++) {
137                  vx = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
138 <                if (vx >= fvxr-1) vx--;
138 >                while (vx >= fvxr-1) vx--;
139                  dx -= (double)vx;
140                  for (i = 0; i < 3; i++) {
141                          lv = (1.-dy)*colval(veilscan(vy)[vx],i) +
# Line 137 | Line 146 | int    y;
146                                          (1.-dx)*lv + dx*uv;
147                  }
148          }
149 + }
150 +
151 +
152 + /****************** ACUITY STUFF *******************/
153 +
154 + typedef struct {
155 +        short   sampe;          /* sample area size (exponent of 2) */
156 +        short   nscans;         /* number of scanlines in this bar */
157 +        int     len;            /* individual scanline length */
158 +        int     nread;          /* number of scanlines loaded */
159 +        COLOR   *sdata;         /* scanbar data */
160 + } SCANBAR;
161 +
162 + #define bscan(sb,y)     ((COLOR *)(sb)->sdata+((y)%(sb)->nscans)*(sb)->len)
163 +
164 + SCANBAR *rootbar;               /* root scan bar (lowest resolution) */
165 +
166 + float   *inpacuD;               /* input acuity data (cycles/degree) */
167 +
168 + #define tsampr(x,y)     inpacuD[(y)*fvxr+(x)]
169 +
170 +
171 + double
172 + hacuity(La)                     /* return visual acuity in cycles/degree */
173 + double  La;
174 + {
175 +                                        /* functional fit */
176 +        return(17.25*atan(1.4*log10(La) + 0.35) + 25.72);
177 + }
178 +
179 +
180 + COLOR *
181 + getascan(sb, y)                 /* find/read scanline y for scanbar sb */
182 + register SCANBAR        *sb;
183 + int     y;
184 + {
185 +        register COLOR  *sl0, *sl1, *mysl;
186 +        register int    i;
187 +
188 +        if (y < sb->nread - sb->nscans)                 /* too far back? */
189 +                return(NULL);
190 +        for ( ; y >= sb->nread; sb->nread++) {          /* read as necessary */
191 +                mysl = bscan(sb, sb->nread);
192 +                if (sb->sampe == 0) {
193 +                        if (freadscan(mysl, sb->len, infp) < 0) {
194 +                                fprintf(stderr, "%s: %s: scanline read error\n",
195 +                                                progname, infn);
196 +                                exit(1);
197 +                        }
198 +                } else {
199 +                        sl0 = getascan(sb+1, 2*y);
200 +                        if (sl0 == NULL)
201 +                                return(NULL);
202 +                        sl1 = getascan(sb+1, 2*y+1);
203 +                        for (i = 0; i < sb->len; i++) {
204 +                                copycolor(mysl[i], sl0[2*i]);
205 +                                addcolor(mysl[i], sl0[2*i+1]);
206 +                                addcolor(mysl[i], sl1[2*i]);
207 +                                addcolor(mysl[i], sl1[2*i+1]);
208 +                                scalecolor(mysl[i], 0.25);
209 +                        }
210 +                }
211 +        }
212 +        return(bscan(sb, y));
213 + }
214 +
215 +
216 + acuscan(scln, y)                /* get acuity-sampled scanline */
217 + COLOR   *scln;
218 + int     y;
219 + {
220 +        double  sr;
221 +        double  dx, dy;
222 +        int     ix, iy;
223 +        register int    x;
224 +                                        /* compute foveal y position */
225 +        iy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
226 +        while (iy >= fvyr-1) iy--;
227 +        dy -= (double)iy;
228 +        for (x = 0; x < scanlen(&inpres); x++) {
229 +                                        /* compute foveal x position */
230 +                ix = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
231 +                while (ix >= fvxr-1) ix--;
232 +                dx -= (double)ix;
233 +                                        /* interpolate sample rate */
234 +                sr = (1.-dy)*((1.-dx)*tsampr(ix,iy) + dx*tsampr(ix+1,iy)) +
235 +                        dy*((1.-dx)*tsampr(ix,iy+1) + dx*tsampr(ix+1,iy+1));
236 +
237 +                acusample(scln[x], x, y, sr);   /* compute sample */
238 +        }
239 + }
240 +
241 +
242 + acusample(col, x, y, sr)        /* interpolate sample at (x,y) using rate sr */
243 + COLOR   col;
244 + int     x, y;
245 + double  sr;
246 + {
247 +        COLOR   c1;
248 +        double  d;
249 +        register SCANBAR        *sb0;
250 +
251 +        for (sb0 = rootbar; sb0->sampe != 0 && 1<<sb0[1].sampe > sr; sb0++)
252 +                ;
253 +        ascanval(col, x, y, sb0);
254 +        if (sb0->sampe == 0)            /* don't extrapolate highest */
255 +                return;
256 +        ascanval(c1, x, y, sb0+1);
257 +        d = ((1<<sb0->sampe) - sr)/(1<<sb0[1].sampe);
258 +        scalecolor(col, 1.-d);
259 +        scalecolor(c1, d);
260 +        addcolor(col, c1);
261 + }
262 +
263 +
264 + ascanval(col, x, y, sb)         /* interpolate scanbar at orig. coords (x,y) */
265 + COLOR   col;
266 + int     x, y;
267 + SCANBAR *sb;
268 + {
269 +        COLOR   *sl0, *sl1, c1, c1y;
270 +        double  dx, dy;
271 +        int     ix, iy;
272 +
273 +        if (sb->sampe == 0) {           /* no need to interpolate */
274 +                sl0 = getascan(sb, y);
275 +                copycolor(col, sl0[x]);
276 +                return;
277 +        }
278 +                                        /* compute coordinates for sb */
279 +        ix = dx = (x+.5)/(1<<sb->sampe) - .5;
280 +        while (ix >= sb->len-1) ix--;
281 +        dx -= (double)ix;
282 +        iy = dy = (y+.5)/(1<<sb->sampe) - .5;
283 +        while (iy >= (numscans(&inpres)>>sb->sampe)-1) iy--;
284 +        dy -= (double)iy;
285 +                                        /* get scanlines */
286 +        sl0 = getascan(sb, iy);
287 + #ifdef DEBUG
288 +        if (sl0 == NULL)
289 +                error(INTERNAL, "cannot backspace in ascanval");
290 + #endif
291 +        sl1 = getascan(sb, iy+1);
292 +                                        /* 2D linear interpolation */
293 +        copycolor(col, sl0[ix]);
294 +        scalecolor(col, 1.-dx);
295 +        copycolor(c1, sl0[ix+1]);
296 +        scalecolor(c1, dx);
297 +        addcolor(col, c1);
298 +        copycolor(c1y, sl1[ix]);
299 +        scalecolor(c1y, 1.-dx);
300 +        copycolor(c1, sl1[ix+1]);
301 +        scalecolor(c1, dx);
302 +        addcolor(c1y, c1);
303 +        scalecolor(col, 1.-dy);
304 +        scalecolor(c1y, dy);
305 +        addcolor(col, c1y);
306 +        for (ix = 0; ix < 3; ix++)      /* make sure no negative */
307 +                if (colval(col,ix) < 0.)
308 +                        colval(col,ix) = 0.;
309 + }
310 +
311 +
312 + SCANBAR *
313 + sballoc(se, ns, sl)             /* allocate scanbar */
314 + int     se;             /* sampling rate exponent */
315 + int     ns;             /* number of scanlines */
316 + int     sl;             /* original scanline length */
317 + {
318 +        SCANBAR *sbarr;
319 +        register SCANBAR        *sb;
320 +
321 +        sbarr = sb = (SCANBAR *)malloc((se+1)*sizeof(SCANBAR));
322 +        if (sb == NULL)
323 +                syserror("malloc");
324 +        do {
325 +                sb->sampe = se;
326 +                sb->len = sl>>se;
327 +                sb->nscans = ns;
328 +                sb->sdata = (COLOR *)malloc(sb->len*ns*sizeof(COLOR));
329 +                if (sb->sdata == NULL)
330 +                        syserror("malloc");
331 +                sb->nread = 0;
332 +                ns <<= 1;
333 +                sb++;
334 +        } while (--se >= 0);
335 +        return(sbarr);
336 + }
337 +
338 +
339 + initacuity()                    /* initialize variable acuity sampling */
340 + {
341 +        FVECT   diffx, diffy, cp;
342 +        double  omega, maxsr;
343 +        register int    x, y, i;
344 +
345 +        compraydir();                   /* compute ray directions */
346 +
347 +        inpacuD = (float *)malloc(fvxr*fvyr*sizeof(float));
348 +        if (inpacuD == NULL)
349 +                syserror("malloc");
350 +        maxsr = 1.;                     /* compute internal sample rates */
351 +        for (y = 1; y < fvyr-1; y++)
352 +                for (x = 1; x < fvxr-1; x++) {
353 +                        for (i = 0; i < 3; i++) {
354 +                                diffx[i] = 0.5*fvxr/scanlen(&inpres) *
355 +                                                (rdirscan(y)[x+1][i] -
356 +                                                rdirscan(y)[x-1][i]);
357 +                                diffy[i] = 0.5*fvyr/numscans(&inpres) *
358 +                                                (rdirscan(y+1)[x][i] -
359 +                                                rdirscan(y-1)[x][i]);
360 +                        }
361 +                        fcross(cp, diffx, diffy);
362 +                        omega = 0.5 * sqrt(DOT(cp,cp));
363 +                        if (omega <= FTINY*FTINY)
364 +                                tsampr(x,y) = 1.;
365 +                        else if ((tsampr(x,y) = PI/180. / sqrt(omega) /
366 +                                        hacuity(plum(fovscan(y)[x]))) > maxsr)
367 +                                maxsr = tsampr(x,y);
368 +                }
369 +                                        /* copy perimeter (easier) */
370 +        for (x = 1; x < fvxr-1; x++) {
371 +                tsampr(x,0) = tsampr(x,1);
372 +                tsampr(x,fvyr-1) = tsampr(x,fvyr-2);
373 +        }
374 +        for (y = 0; y < fvyr; y++) {
375 +                tsampr(0,y) = tsampr(1,y);
376 +                tsampr(fvxr-1,y) = tsampr(fvxr-2,y);
377 +        }
378 +                                        /* initialize with next power of two */
379 +        rootbar = sballoc((int)(log(maxsr)/log(2.))+1, 2, scanlen(&inpres));
380   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines