ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmglare.c
Revision: 2.1
Committed: Mon Sep 9 17:19:51 2019 UTC (3 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R3, HEAD
Log Message:
Added Nathaniel Jones' dcglare utility

File Contents

# Content
1 /*
2 * cmglare.c - routines for calculating glare autonomy.
3 *
4 * N. Jones
5 */
6
7 /*
8 * Copyright (c) 2017-2019 Nathaniel Jones
9 *
10 * Permission is hereby granted, free of charge, to any person obtaining a copy
11 * of this software and associated documentation files (the "Software"), to deal
12 * in the Software without restriction, including without limitation the rights
13 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
14 * copies of the Software, and to permit persons to whom the Software is
15 * furnished to do so, subject to the following conditions:
16 *
17 * The above copyright notice and this permission notice shall be included in all
18 * copies or substantial portions of the Software.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
21 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
23 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
25 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
26 * SOFTWARE.
27 */
28
29 #include "rtmath.h"
30 #include "cmglare.h"
31
32 #define LUMINOUS_EFFICACY 179 /* lumens per Watt */
33 #define LUMINANCE_THRESHOLD 100 /* minimum threshold that will be interpreted as luminance rather than ratio between Ev and glare source */
34
35 #define ANGLE(u,v) acos(DOT((u),(v))) // TODO need to normalize?
36
37 typedef struct reinhart_sky {
38 //int mf; /* Linear divisions per Tregenza patch. */
39 int rings; /* Number of rings of sky patches. */
40 //int patches; /* Number of sky patches. */
41 double ringElevationAngle; /* Angle in radians between rings. */
42 int *patchesPerRow; /* Sky patches per row. */
43 int *firstPatchIndex; /* Index of first sky patch in each row. */
44 double *solidAngle; /* Solid angle of each patch in each row. */
45 } ReinhartSky;
46
47 static const int tnaz[] = { 30, 30, 24, 24, 18, 12, 6 }; /* Number of patches per row */
48
49 extern char* progname;
50
51 static ReinhartSky* make_sky(const CMATRIX *smx)
52 {
53 int mf, count, i, j;
54 double remaining = PI / 2;
55
56 /* Check sky partitionss */
57 switch (smx->nrows) { // nrows is Reinhart sky subdivisoins plus one for miss
58 case 2: mf = 0; break;
59 case 146: mf = 1; break;
60 case 578: mf = 2; break;
61 case 1298: mf = 3; break;
62 case 2306: mf = 4; break;
63 case 3602: mf = 5; break;
64 case 5186: mf = 6; break;
65 case 7058: mf = 7; break;
66 case 9218: mf = 8; break;
67 case 11666: mf = 9; break;
68 case 14402: mf = 10; break;
69 case 17426: mf = 11; break;
70 case 20738: mf = 12; break;
71 default:
72 fprintf(stderr,
73 "%s: unknown number of sky patches %d\n",
74 progname, smx->nrows);
75 return NULL;
76 }
77
78 /* Allocate sky */
79 ReinhartSky *sky = (ReinhartSky*)malloc(sizeof(ReinhartSky));
80 if (!sky) goto memerr;
81
82 /* Calculate patches per row */
83 sky->rings = 7 * mf + 1;
84 sky->patchesPerRow = (int*)malloc(sky->rings * sizeof(int));
85 sky->firstPatchIndex = (int*)malloc(sky->rings * sizeof(int));
86 if (!sky->patchesPerRow || !sky->firstPatchIndex) goto memerr;
87 count = 1; // The below horizon patch
88 for (i = 0; i < 7; i++)
89 for (j = 0; j < mf; j++) {
90 sky->firstPatchIndex[i * mf + j] = count;
91 count += sky->patchesPerRow[i * mf + j] = tnaz[i] * mf;
92 }
93 sky->firstPatchIndex[7 * mf] = count;
94 sky->patchesPerRow[7 * mf] = 1;
95 //sky->patches = count + 1;
96
97 /* Calculate solid angle of patches in each row */
98 sky->solidAngle = (double*)malloc(sky->rings * sizeof(double));
99 if (!sky->solidAngle) goto memerr;
100 sky->ringElevationAngle = PI / (2 * sky->rings - 1);
101 sky->solidAngle[0] = 2 * PI;
102 for (i = 1; i < sky->rings; i++) {
103 remaining -= sky->ringElevationAngle;
104 sky->solidAngle[i] = 2 * PI * (1 - cos(remaining)); // solid angle of cap
105 sky->solidAngle[i - 1] -= sky->solidAngle[i];
106 sky->solidAngle[i - 1] /= sky->patchesPerRow[i - 1];
107 }
108
109 return sky;
110
111 memerr:
112 fprintf(stderr,
113 "%s: out of memory\n",
114 progname);
115 return NULL;
116 }
117
118 static void free_sky(ReinhartSky *sky)
119 {
120 if (sky) {
121 free(sky->patchesPerRow);
122 free(sky->firstPatchIndex);
123 free(sky->solidAngle);
124 free(sky);
125 }
126 }
127
128 static void get_patch_direction(const ReinhartSky *sky, const int patch, FVECT target)
129 {
130 //if (patch >= sky->patches || patch < 0) throw new RuntimeException("Illegal patch " + patch);
131 if (!patch) {
132 target[0] = target[1] = 0;
133 target[2] = -1;
134 return; // Ignore below horizon?
135 }
136 int row = sky->rings - 1;
137 while (patch < sky->firstPatchIndex[row]) row--;
138 const double alt = PI / 2 - sky->ringElevationAngle * (row - (sky->rings - 1));
139 const double azi = 2 * PI * (patch - sky->firstPatchIndex[row]) / sky->patchesPerRow[row];
140 const double cos_alt = cos(alt);
141 target[0] = cos_alt * -sin(azi);
142 target[1] = cos_alt * -cos(azi);
143 target[2] = sin(alt);
144 }
145
146 static double get_patch_solid_angle(const ReinhartSky *sky, const int patch, const double cos_theta)
147 {
148 //if (patch >= sky->patches || patch < 0) throw new RuntimeException("Illegal patch " + patch);
149 if (!patch) return 2 * (PI - 2 * acos(cos_theta)); // Solid angle overlap between visible hemisphere and ground hemisphere
150 int row = sky->rings - 1;
151 while (patch < sky->firstPatchIndex[row]) row--;
152 return sky->solidAngle[row];
153 }
154
155 static double get_guth(const FVECT dir, const FVECT forward, const FVECT up)
156 {
157 double posindex;
158 FVECT hv, temp;
159 VCROSS(hv, forward, up);
160 normalize(hv);
161 VCROSS(temp, forward, hv);
162 double phi = ANGLE(dir, temp) - PI / 2;
163
164 /* Guth model, equation from IES lighting handbook */
165 if (phi >= 0) {
166 double sigma = ANGLE(dir, forward);
167 VSUM(hv, forward, dir, 1 / DOT(dir, forward));
168 normalize(hv);
169 double tau = ANGLE(hv, up);
170 tau *= 180.0 / PI;
171 sigma *= 180.0 / PI;
172
173 if (sigma <= 0)
174 sigma = -sigma;
175
176 posindex = exp((35.2 - 0.31889 * tau - 1.22 * exp(-2.0 * tau / 9.0)) / 1000.0 * sigma + (21.0 + 0.26667 * tau - 0.002963 * tau * tau) / 100000.0 * sigma * sigma);
177 }
178 /* below line of sight, using Iwata model */
179 else {
180 double teta = PI / 2 - ANGLE(dir, hv);
181
182 if (teta == 0.0)
183 teta = FTINY;
184
185 double fact = 0.8;
186 double d = 1 / tan(phi);
187 double s = tan(teta) / tan(phi);
188 double r = sqrt((1 + s * s) / (d * d));
189 if (r > 0.6)
190 fact = 1.2;
191 if (r > 3)
192 r = 3.0;
193
194 posindex = 1.0 + fact * r;
195 }
196 if (posindex > 16)
197 posindex = 16.0;
198
199 return posindex;
200 }
201
202 float* cm_glare(const CMATRIX *dcmx, const CMATRIX *evmx, const CMATRIX *smx, const int *occupied, const double dgp_limit, const double dgp_threshold, const FVECT *views, const FVECT dir, const FVECT up)
203 {
204 int p, t, c;
205 int hourly_output = dgp_limit < 0;
206 float *dgp_list;
207 ReinhartSky *sky;
208 FVECT vdir;
209
210 /* Check consistancy */
211 if ((dcmx->nrows != evmx->nrows) | (dcmx->ncols != smx->nrows) | (evmx->ncols != smx->ncols)) {
212 fprintf(stderr,
213 "%s: inconsistant matrix dimensions: dc(%d, %d) ev(%d, %d) s(%d, %d)\n",
214 progname, dcmx->nrows, dcmx->ncols, evmx->nrows, evmx->ncols, smx->nrows, smx->ncols);
215 return NULL;
216 }
217
218 /* Create output buffer */
219 dgp_list = (float*)malloc(evmx->nrows * (hourly_output ? evmx->ncols : 1) * sizeof(float));
220 if (!dgp_list) {
221 fprintf(stderr,
222 "%s: out of memory in cm_glare()\n",
223 progname);
224 return NULL;
225 }
226
227 /* Create sky */
228 sky = make_sky(smx);
229 if (!sky) return NULL;
230
231 /* Calculate glare limit */
232 double ev_max = -1;
233 if (!hourly_output) {
234 ev_max = (dgp_limit - 0.16) / 5.87e-5;
235 if (ev_max < 0) ev_max = 0;
236 }
237
238 /* For each position and direction */
239 if (!views) VCOPY(vdir, dir);
240 for (p = 0; p < evmx->nrows; p++) {
241 /* For each time step */
242 int occupied_hours = 0;
243 int glare_hours = 0;
244 if (views) VCOPY(vdir, views[p]);
245 for (t = 0; t < evmx->ncols; t++) {
246 if (!occupied[t]) {
247 /* Not occupied */
248 if (hourly_output) dgp_list[p * evmx->ncols + t] = 0.0f;
249 }
250 else {
251 /* Occupied */
252 double illum = LUMINOUS_EFFICACY * bright(cm_lval(evmx, p, t));
253 occupied_hours++;
254
255 if (illum <= FTINY) {
256 /* No light, so no glare */
257 if (hourly_output) dgp_list[p * evmx->ncols + t] = 0.0f;
258 }
259 else if ((illum >= ev_max) & (!hourly_output)) {
260 /* Guarangeed glare */
261 glare_hours++;
262 }
263 else {
264 /* Calculate enhanced simplified daylight glare probability */
265 double sum = 0.0;
266 FVECT patch_normal;
267 for (c = 0; c < smx->nrows; c++) {
268 const double dc = bright(cm_lval(dcmx, p, c));
269 if (dc > 0) {
270 get_patch_direction(sky, c, patch_normal);
271 if (!c) {
272 /* Direction toward the center of the visible ground */
273 VADD(patch_normal, patch_normal, vdir);
274 if (normalize(patch_normal) == 0) continue;
275 }
276 const double cos_theta = DOT(vdir, patch_normal);
277 if (cos_theta <= FTINY) continue;
278 const double s = LUMINOUS_EFFICACY * bright(cm_lval(smx, c, t));
279 const double omega = get_patch_solid_angle(sky, c, cos_theta);
280 const double patch_luminance = (dc * s) / (omega * cos_theta);
281
282 double min_patch_luminance = dgp_threshold;
283 if (dgp_threshold < LUMINANCE_THRESHOLD)
284 min_patch_luminance *= illum / PI; // TODO should use average luminance, not illuminance
285 if (patch_luminance < min_patch_luminance) continue;
286 const double P = get_guth(patch_normal, vdir, up);
287 sum += (patch_luminance * patch_luminance * omega) / (P * P);
288 }
289 }
290
291 //double dgp = 5.87e-5 * illum + 0.092 * log10(1 + dgp / pow(illum, 1.87)) + 0.159;
292 double eDGPs = 5.87e-5 * illum + 0.0918 * log10(1 + sum / pow(illum, 1.87)) + 0.16;
293 if (illum < 1000) /* low light correction */
294 eDGPs *= exp(0.024 * illum - 4) / (1 + exp(0.024 * illum - 4));
295 //eDGPs /= 1.1 - 0.5 * age / 100.0; /* age correction */
296 if (eDGPs > 1.0) eDGPs = 1.0;
297
298 if (hourly_output)
299 dgp_list[p * evmx->ncols + t] = (float)eDGPs;
300 else if (eDGPs >= dgp_limit)
301 glare_hours++;
302 }
303 }
304 }
305 if (!hourly_output) {
306 /* Save glare autonomy */
307 dgp_list[p] = (float)(occupied_hours - glare_hours) / occupied_hours;
308 }
309 }
310
311 free_sky(sky);
312
313 return dgp_list;
314 }
315
316 static int getvec(FVECT vec, const int dtype, FILE *fp) /* get a vector from fp */
317 {
318 static float vf[3];
319 static double vd[3];
320 char buf[32];
321 int i;
322
323 switch (dtype) {
324 case DTascii:
325 for (i = 0; i < 3; i++) {
326 if (fgetword(buf, sizeof(buf), fp) == NULL ||
327 !isflt(buf))
328 return(-1);
329 vec[i] = atof(buf);
330 }
331 break;
332 case DTfloat:
333 if (getbinary(vf, sizeof(float), 3, fp) != 3)
334 return(-1);
335 VCOPY(vec, vf);
336 break;
337 case DTdouble:
338 if (getbinary(vd, sizeof(double), 3, fp) != 3)
339 return(-1);
340 VCOPY(vec, vd);
341 break;
342 default:
343 fprintf(stderr,
344 "%s: botched input format\n",
345 progname);
346 return(-1);
347 }
348 return(0);
349 }
350
351 int cm_load_schedule(const int count, int* schedule, FILE *fp)
352 {
353 char buf[512];
354 char *cp;
355 char *comma;
356 double val;
357 int i = 0;
358
359 while (fgetline(buf, sizeof(buf), fp) != NULL) {
360 if (buf[0] == '#') continue; // Comment line
361 comma = NULL;
362 for (cp = buf; *cp; cp++) {
363 /* If there are multiple commas, assume the value is after the last comma */
364 if (*cp == ',') {
365 comma = cp; /* Record position of last comma */
366 }
367 }
368 if (comma)
369 val = atof(comma + 1);
370 else
371 val = atof(buf);
372
373 if (i < count) {
374 /* Add the value to the schedule */
375 schedule[i++] = (val > 0);
376 }
377 else {
378 fprintf(stderr,
379 "%s: too many schedule entries\n",
380 progname);
381 return(1);
382 }
383 }
384 fclose(fp);
385
386 if (i < count) {
387 fprintf(stderr,
388 "%s: too few schedule entries\n",
389 progname);
390 return(-1);
391 }
392
393 return 0;
394 }
395
396 FVECT* cm_load_views(const int nrows, const int dtype, FILE *fp)
397 {
398 int i;
399 double d;
400 FVECT orig;
401
402 FVECT *views = (FVECT*)malloc(nrows * sizeof(FVECT));
403 if (!views) {
404 fprintf(stderr,
405 "%s: out of memory in cm_load_views()\n",
406 progname);
407 return NULL;
408 }
409
410 for (i = 0; i < nrows; i++) {
411 if (getvec(orig, dtype, fp) | getvec(views[i], dtype, fp)) {
412 fprintf(stderr,
413 "%s: unexpected end of input, missing %d entries\n",
414 progname, i);
415 return NULL;
416 }
417
418 d = normalize(views[i]);
419 if (d == 0.0) { /* zero ==> flush */
420 fprintf(stderr,
421 "%s: zero length direction detected\n",
422 progname);
423 return NULL;
424 }
425 }
426
427 return views;
428 }
429
430 int cm_write_glare(const float *mp, const int nrows, const int ncols, const int dtype, FILE *fp)
431 {
432 static const char tabEOL[2] = { '\t', '\n' };
433 int r, c;
434 double dc[1];
435
436 switch (dtype) {
437 case DTascii:
438 for (r = 0; r < nrows; r++)
439 for (c = 0; c < ncols; c++, mp++)
440 fprintf(fp, "%.6e%c",
441 mp[0],
442 tabEOL[c >= ncols - 1]);
443 break;
444 case DTfloat:
445 r = ncols*nrows;
446 while (r > 0) {
447 c = putbinary(mp, sizeof(float), r, fp);
448 if (c <= 0)
449 return(0);
450 mp += c;
451 r -= c;
452 }
453 break;
454 case DTdouble:
455 r = ncols*nrows;
456 while (r--) {
457 dc[0] = mp[0];
458 if (putbinary(dc, sizeof(double), 1, fp) != 1)
459 return(0);
460 mp++;
461 }
462 break;
463 default:
464 fputs("Unsupported data type in cm_write_glare()!\n", stderr);
465 return(0);
466 }
467 return(fflush(fp) == 0);
468 }