ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/px/pcond4.c

Comparing ray/src/px/pcond4.c (file contents):
Revision 3.9 by greg, Thu Oct 10 16:36:39 1996 UTC vs.
Revision 3.16 by greg, Sat Feb 22 02:07:27 2003 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1996 Regents of the University of California */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ LBL";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		* Routines for veiling glare and loss of acuity.
6		*/
#	Line 14 | Line 11 | static char SCCSid[] = "$SunId$ LBL";
11
12		#define VADAPT          0.08            /* fraction of adaptation from veil */
13
14	<	extern COLOR    *fovimg;                /* foveal (1 degree) averaged image */
18	<	extern short    fvxr, fvyr;             /* foveal image resolution */
14	>	static COLOR    *veilimg = NULL;        /* veiling image */
15
20	–	#define fovscan(y)      (fovimg+(y)*fvxr)
21	–
22	–	static COLOR    *veilimg;               /* veiling image */
23	–
16		#define veilscan(y)     (veilimg+(y)*fvxr)
17
18		static float    (*raydir)[3] = NULL;    /* ray direction for each pixel */
#	Line 41 | Line 33 | compraydir()                           /* compute ray directions */
33		syserror("malloc");
34
35		for (y = 0; y < fvyr; y++) {
36	<	switch (inpres.or) {
36	>	switch (inpres.rt) {
37		case YMAJOR: case YMAJOR|XDECR:
38		v = (y+.5)/fvyr; break;
39		case YMAJOR|YDECR: case YMAJOR|YDECR|XDECR:
#	Line 52 | Line 44 | compraydir()                           /* compute ray directions */
44		h = 1. - (y+.5)/fvyr; break;
45		}
46		for (x = 0; x < fvxr; x++) {
47	<	switch (inpres.or) {
47	>	switch (inpres.rt) {
48		case YMAJOR: case YMAJOR|YDECR:
49		h = (x+.5)/fvxr; break;
50		case YMAJOR|XDECR: case YMAJOR|XDECR|YDECR:
#	Line 83 | Line 75 | compveil()                             /* compute veiling image */
75		COLOR   ctmp, vsum;
76		int     px, py;
77		register int    x, y;
78	+
79	+	if (veilimg != NULL)            /* already done? */
80	+	return;
81		/* compute ray directions */
82		compraydir();
83		/* compute veil image */
#	Line 100 | Line 95 | compveil()                             /* compute veiling image */
95		rdirscan(y)[x]);
96		if (t2 <= FTINY) continue;
97		/*      use approximation instead
98	<	t2 = acos(t2);
99	<	t2 = 1./(t2*t2);
98	>	t3 = acos(t2);
99	>	t2 = t2/(t3*t3);
100		*/
101	<	t2 = .5 / (1. - t2);
101	>	t2 *= .5 / (1. - t2);
102		copycolor(ctmp, fovscan(y)[x]);
103		scalecolor(ctmp, t2);
104		addcolor(vsum, ctmp);
105		t2sum += t2;
106		}
107		/* VADAPT of original is subtracted in addveil() */
108	<	scalecolor(vsum, VADAPT/t2sum);
108	>	if (t2sum > FTINY)
109	>	scalecolor(vsum, VADAPT/t2sum);
110		copycolor(veilscan(py)[px], vsum);
111		}
112	+	/* modify FOV sample image */
113	+	for (y = 0; y < fvyr; y++)
114	+	for (x = 0; x < fvxr; x++) {
115	+	scalecolor(fovscan(y)[x], 1.-VADAPT);
116	+	addcolor(fovscan(y)[x], veilscan(y)[x]);
117	+	}
118	+	comphist();                     /* recompute histogram */
119		}
120
121
122	+	#if ADJ_VEIL
123	+	/*
124	+	* The following veil adjustment was added to compensate for
125	+	* the fact that contrast reduction gets confused with veil
126	+	* in the human visual system.  Therefore, we reduce the
127	+	* veil in portions of the image where our mapping has
128	+	* already reduced contrast below the target value.
129	+	* This gets called after the intial veil has been computed
130	+	* and added to the foveal image, and the mapping has been
131	+	* determined.
132	+	*/
133	+	adjveil()                               /* adjust veil image */
134	+	{
135	+	float   *crfptr = crfimg;
136	+	COLOR   *fovptr = fovimg;
137	+	COLOR   *veilptr = veilimg;
138	+	double  s2nits = 1./inpexp;
139	+	double  vl, vl2, fovl, vlsum;
140	+	double  deltavc[3];
141	+	int     i, j;
142	+
143	+	if (lumf == rgblum)
144	+	s2nits *= WHTEFFICACY;
145	+
146	+	for (i = fvxr*fvyr; i--; crfptr++, fovptr++, veilptr++) {
147	+	if (crfptr[0] >= 0.95)
148	+	continue;
149	+	vl = plum(veilptr[0]);
150	+	fovl = (plum(fovptr[0]) - vl) * (1./(1.-VADAPT));
151	+	if (vl <= 0.05*fovl)
152	+	continue;
153	+	vlsum = vl;
154	+	for (j = 2; j < 11; j++) {
155	+	vlsum += crfptr[0]*vl - (1.0 - crfptr[0])*fovl;
156	+	vl2 = vlsum / (double)j;
157	+	if (vl2 < 0.0)
158	+	vl2 = 0.0;
159	+	crfptr[0] = crfactor(fovl + vl2);
160	+	}
161	+	/* desaturation code causes color fringes at this level */
162	+	for (j = 3; j--; ) {
163	+	double  vc = colval(veilptr[0],j);
164	+	double  fovc = (colval(fovptr[0],j) - vc) *
165	+	(1./(1.-VADAPT));
166	+	deltavc[j] = (1.-crfptr[0])*(fovl/s2nits - fovc);
167	+	if (vc + deltavc[j] < 0.0)
168	+	break;
169	+	}
170	+	if (j < 0)
171	+	addcolor(veilptr[0], deltavc);
172	+	else
173	+	scalecolor(veilptr[0], vl2/vl);
174	+	}
175	+	smoothveil();                   /* smooth our result */
176	+	}
177	+
178	+
179	+	smoothveil()                            /* smooth veil image */
180	+	{
181	+	COLOR   *nveilimg;
182	+	COLOR   *ovptr, *nvptr;
183	+	int     x, y, i;
184	+
185	+	nveilimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR));
186	+	if (nveilimg == NULL)
187	+	return;
188	+	for (y = 1; y < fvyr-1; y++) {
189	+	ovptr = veilimg + y*fvxr + 1;
190	+	nvptr = nveilimg + y*fvxr + 1;
191	+	for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
192	+	for (i = 3; i--; )
193	+	nvptr[0][i] = 0.5 * ovptr[0][i]
194	+	+ (1./12.) *
195	+	(ovptr[-1][i] + ovptr[-fvxr][i] +
196	+	ovptr[1][i] + ovptr[fvxr][i])
197	+	+ (1./24.) *
198	+	(ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i] +
199	+	ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
200	+	}
201	+	ovptr = veilimg + 1;
202	+	nvptr = nveilimg + 1;
203	+	for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
204	+	for (i = 3; i--; )
205	+	nvptr[0][i] = 0.5 * ovptr[0][i]
206	+	+ (1./9.) *
207	+	(ovptr[-1][i] + ovptr[1][i] + ovptr[fvxr][i])
208	+	+ (1./12.) *
209	+	(ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
210	+	ovptr = veilimg + (fvyr-1)*fvxr + 1;
211	+	nvptr = nveilimg + (fvyr-1)*fvxr + 1;
212	+	for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
213	+	for (i = 3; i--; )
214	+	nvptr[0][i] = 0.5 * ovptr[0][i]
215	+	+ (1./9.) *
216	+	(ovptr[-1][i] + ovptr[1][i] + ovptr[-fvxr][i])
217	+	+ (1./12.) *
218	+	(ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i]);
219	+	ovptr = veilimg + fvxr;
220	+	nvptr = nveilimg + fvxr;
221	+	for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
222	+	for (i = 3; i--; )
223	+	nvptr[0][i] = 0.5 * ovptr[0][i]
224	+	+ (1./9.) *
225	+	(ovptr[-fvxr][i] + ovptr[1][i] + ovptr[fvxr][i])
226	+	+ (1./12.) *
227	+	(ovptr[-fvxr+1][i] + ovptr[fvxr+1][i]);
228	+	ovptr = veilimg + fvxr - 1;
229	+	nvptr = nveilimg + fvxr - 1;
230	+	for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
231	+	for (i = 3; i--; )
232	+	nvptr[0][i] = 0.5 * ovptr[0][i]
233	+	+ (1./9.) *
234	+	(ovptr[-fvxr][i] + ovptr[-1][i] + ovptr[fvxr][i])
235	+	+ (1./12.) *
236	+	(ovptr[-fvxr-1][i] + ovptr[fvxr-1][i]);
237	+	for (i = 3; i--; ) {
238	+	nveilimg[0][i] = veilimg[0][i];
239	+	nveilimg[fvxr-1][i] = veilimg[fvxr-1][i];
240	+	nveilimg[(fvyr-1)*fvxr][i] = veilimg[(fvyr-1)*fvxr][i];
241	+	nveilimg[fvyr*fvxr-1][i] = veilimg[fvyr*fvxr-1][i];
242	+	}
243	+	free((void *)veilimg);
244	+	veilimg = nveilimg;
245	+	}
246	+	#endif
247	+
248		addveil(sl, y)                          /* add veil to scanline */
249		COLOR   *sl;
250		int     y;
#	Line 166 | Line 295 | float  *inpacuD;               /* input acuity data (cycles/degree)
295		double
296		hacuity(La)                     /* return visual acuity in cycles/degree */
297		double  La;
298	<	{                               /* data due to S. Shaler (we should fit it!) */
299	<	#define NPOINTS 20
300	<	static float    l10lum[NPOINTS] = {
172	<	-3.10503,-2.66403,-2.37703,-2.09303,-1.64403,-1.35803,
173	<	-1.07403,-0.67203,-0.38503,-0.10103,0.29397,0.58097,0.86497,
174	<	1.25697,1.54397,1.82797,2.27597,2.56297,2.84697,3.24897
175	<	};
176	<	static float    resfreq[NPOINTS] = {
177	<	2.09,3.28,3.79,4.39,6.11,8.83,10.94,18.66,23.88,31.05,37.42,
178	<	37.68,41.60,43.16,45.30,47.00,48.43,48.32,51.06,51.09
179	<	};
180	<	double  l10La;
181	<	register int    i;
182	<	/* check limits */
183	<	if (La <= 7.85e-4)
184	<	return(resfreq[0]);
185	<	if (La >= 1.78e3)
186	<	return(resfreq[NPOINTS-1]);
187	<	/* interpolate data */
188	<	l10La = log10(La);
189	<	for (i = 0; i < NPOINTS-2 && l10lum[i+1] <= l10La; i++)
190	<	;
191	<	return( ( (l10lum[i+1] - l10La)*resfreq[i] +
192	<	(l10La - l10lum[i])*resfreq[i+1] ) /
193	<	(l10lum[i+1] - l10lum[i]) );
194	<	#undef NPOINTS
298	>	{
299	>	/* functional fit */
300	>	return(17.25*atan(1.4*log10(La) + 0.35) + 25.72);
301		}
302
303
#	Line 303 | Line 409 | SCANBAR        *sb;
409		/* get scanlines */
410		sl0 = getascan(sb, iy);
411		#ifdef DEBUG
412	<	if (sl0 == NULL) {
413	<	fprintf(stderr, "%s: internal - cannot backspace in ascanval\n",
308	<	progname);
309	<	abort();
310	<	}
412	>	if (sl0 == NULL)
413	>	error(INTERNAL, "cannot backspace in ascanval");
414		#endif
415		sl1 = getascan(sb, iy+1);
416		/* 2D linear interpolation */
#	Line 343 | Line 446 | int    sl;             /* original scanline length */
446		if (sb == NULL)
447		syserror("malloc");
448		do {
346	–	sb->sampe = se;
449		sb->len = sl>>se;
450	+	if (sb->len <= 0)
451	+	continue;
452	+	sb->sampe = se;
453		sb->nscans = ns;
454		sb->sdata = (COLOR *)malloc(sb->len*ns*sizeof(COLOR));
455		if (sb->sdata == NULL)
#	Line 381 | Line 486 | initacuity()                   /* initialize variable acuity sampling
486		}
487		fcross(cp, diffx, diffy);
488		omega = 0.5 * sqrt(DOT(cp,cp));
489	<	if (omega <= FTINY)
489	>	if (omega <= FTINY*FTINY)
490		tsampr(x,y) = 1.;
491		else if ((tsampr(x,y) = PI/180. / sqrt(omega) /
492		hacuity(plum(fovscan(y)[x]))) > maxsr)

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)