ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/genssky.c
Revision: 2.1
Committed: Fri Jul 5 18:04:36 2024 UTC (10 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
feat(genssky): Taoning Wang added utility for generating spectral skies

File Contents

# User Rev Content
1 greg 2.1 // Main function for generating spectral sky
2     // Cloudy sky computed as weight average of clear and cie overcast sky
3    
4     #include "copyright.h"
5     #include "atmos.h"
6     #include "resolu.h"
7     #include "view.h"
8    
9    
10     char *progname;
11    
12     const double ARCTIC_LAT = 67.;
13     const double TROPIC_LAT = 23.;
14     const int SUMMER_START = 4;
15     const int SUMMER_END = 9;
16     const double GNORM = 0.777778;
17    
18     const double D65EFF = 203.; /* standard illuminant D65 */
19    
20     // Mean normalized relative daylight spectra where CCT = 6415K for overcast;
21     const double D6415[NSSAMP] = {0.63231, 1.06171, 1.00779, 1.36423, 1.34133,
22     1.27258, 1.26276, 1.26352, 1.22201, 1.13246,
23     1.0434, 1.05547, 0.98212, 0.94445, 0.9722,
24     0.82387, 0.87853, 0.82559, 0.75111, 0.78925};
25    
26     static inline double wmean2(const double a, const double b, const double x) {
27     return a * (1 - x) + b * x;
28     }
29    
30     static inline double wmean(const double a, const double x, const double b,
31     const double y) {
32     return (a * x + b * y) / (a + b);
33     }
34    
35     static double get_zenith_brightness(const double sundir[3]) {
36     double zenithbr;
37     if (sundir[2] < 0) {
38     zenithbr = 0;
39     } else {
40     zenithbr = (8.6 * sundir[2] + .123) * 1000.0 / D65EFF;
41     }
42     return zenithbr;
43     }
44    
45     // from gensky.c
46     static double get_overcast_brightness(const double dz, const double zenithbr) {
47     double groundbr = zenithbr * GNORM;
48     return wmean(pow(dz + 1.01, 10), zenithbr * (1 + 2 * dz) / 3,
49     pow(dz + 1.01, -10), groundbr);
50     }
51    
52     static void write_rad_file(FILE *fp, const double *sun_radiance,
53     const FVECT sundir, const char skyfile[PATH_MAX],
54     const char grndfile[PATH_MAX]) {
55     if (sundir[2] > 0) {
56     fprintf(fp, "void spectrum sunrad\n0\n0\n22 380 780 ");
57     for (int i = 0; i < NSSAMP; ++i) {
58     fprintf(fp, "%.1f ", sun_radiance[i] * WVLSPAN);
59     }
60     fprintf(fp, "\n\nsunrad light solar\n0\n0\n3 1 1 1\n\n");
61     fprintf(fp, "solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0],
62     sundir[1], sundir[2]);
63     }
64     fprintf(fp,
65     "void specpict skyfunc\n8 noop %s fisheye.cal fish_u fish_v -rx 90 "
66     "-mx\n0\n0\n\n",
67     skyfile);
68     fprintf(fp, "skyfunc glow sky_glow\n0\n0\n4 1 1 1 0\n\n");
69     fprintf(fp, "sky_glow source sky\n0\n0\n4 0 0 1 180\n\n");
70    
71     fprintf(fp,
72     "void specpict grndmap\n8 noop %s fisheye.cal fish_u fish_v -rx -90 "
73     "-my\n0\n0\n\n",
74     grndfile);
75     fprintf(fp, "grndmap glow ground_glow\n0\n0\n4 1 1 1 0\n\n");
76     fprintf(fp, "ground_glow source ground_source\n0\n0\n4 0 0 -1 180\n\n");
77     }
78    
79     static void write_hsr_header(FILE *fp, RESOLU *res) {
80     float wvsplit[4] = {380, 480, 588,
81     780}; // RGB wavelength limits+partitions (nm)
82     newheader("RADIANCE", fp);
83     fputncomp(NSSAMP, fp);
84     fputwlsplit(wvsplit, fp);
85     fputformat(SPECFMT, fp);
86     fputc('\n', fp);
87     fputsresolu(res, fp);
88     }
89    
90     int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_clear,
91     DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
92     const double cloud_cover, const FVECT sundir,
93     const double grefl, const int res, const char *outname) {
94    
95     char radfile[PATH_MAX];
96     char skyfile[PATH_MAX];
97     char grndfile[PATH_MAX];
98     if (!snprintf(radfile, sizeof(radfile), "%s.rad", outname)) {
99     fprintf(stderr, "Error setting rad file name\n");
100     return 0;
101     };
102     if (!snprintf(skyfile, sizeof(skyfile), "%s_sky.hsr", outname)) {
103     fprintf(stderr, "Error setting sky file name\n");
104     return 0;
105     };
106     if (!snprintf(grndfile, sizeof(grndfile), "%s_ground.hsr", outname)) {
107     fprintf(stderr, "Error setting ground file name\n");
108     return 0;
109     }
110     RESOLU rs = {PIXSTANDARD, res, res};
111     FILE *skyfp = fopen(skyfile, "w");
112     FILE *grndfp = fopen(grndfile, "w");
113     write_hsr_header(grndfp, &rs);
114     write_hsr_header(skyfp, &rs);
115     VIEW skyview = {VT_ANG, {0., 0., 0.}, {0., 0., 1.}, {0., 1., 0.}, 1.,
116     180., 180., 0., 0., 0.,
117     0., {0., 0., 0.}, {0., 0., 0.}, 0., 0.};
118     VIEW grndview = {
119     VT_ANG, {0., 0., 0.}, {0., 0., -1.}, {0., 1., 0.}, 1., 180., 180., 0., 0.,
120     0., 0., {0., 0., 0.}, {0., 0., 0.}, 0., 0.};
121     setview(&skyview);
122     setview(&grndview);
123    
124     CNDX[3] = NSSAMP;
125    
126     FVECT view_point = {0, 0, ER};
127     const double radius = VLEN(view_point);
128     const double sun_ct = fdot(view_point, sundir) / radius;
129     for (unsigned int j = 0; j < res; ++j) {
130     for (unsigned int i = 0; i < res; ++i) {
131     RREAL loc[2];
132     FVECT rorg = {0};
133     FVECT rdir_sky = {0};
134     FVECT rdir_grnd = {0};
135     SCOLOR sky_radiance = {0};
136     SCOLOR ground_radiance = {0};
137     SCOLR sky_sclr = {0};
138     SCOLR ground_sclr = {0};
139    
140     pix2loc(loc, &rs, i, j);
141     viewray(rorg, rdir_sky, &skyview, loc[0], loc[1]);
142     viewray(rorg, rdir_grnd, &grndview, loc[0], loc[1]);
143    
144     const double mu_sky = fdot(view_point, rdir_sky) / radius;
145     const double nu_sky = fdot(rdir_sky, sundir);
146    
147     const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
148     const double nu_grnd = fdot(rdir_grnd, sundir);
149    
150     get_sky_radiance(scat_clear, scat1m_clear, radius, mu_sky, sun_ct, nu_sky,
151     sky_radiance);
152     get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
153     view_point, rdir_grnd, radius, mu_grnd, sun_ct,
154     nu_grnd, grefl, sundir, ground_radiance);
155    
156     for (int k = 0; k < NSSAMP; ++k) {
157     sky_radiance[k] *= WVLSPAN;
158     ground_radiance[k] *= WVLSPAN;
159     }
160    
161     if (cloud_cover > 0) {
162     double zenithbr = get_zenith_brightness(sundir);
163     double grndbr = zenithbr * GNORM;
164     double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
165     for (int k = 0; k < NSSAMP; ++k) {
166     sky_radiance[k] =
167     wmean2(sky_radiance[k], skybr * D6415[k], cloud_cover);
168     ground_radiance[k] =
169     wmean2(ground_radiance[k], grndbr * D6415[k], cloud_cover);
170     }
171     }
172    
173     scolor2scolr(sky_sclr, sky_radiance, 20);
174     putbinary(sky_sclr, LSCOLR, 1, skyfp);
175    
176     scolor2scolr(ground_sclr, ground_radiance, 20);
177     putbinary(ground_sclr, LSCOLR, 1, grndfp);
178     }
179     }
180     fclose(skyfp);
181     fclose(grndfp);
182    
183     // Get solar radiance
184     double sun_radiance[NSSAMP] = {0};
185     get_solar_radiance(tau_clear, scat_clear, scat1m_clear, sundir, radius,
186     sun_ct, sun_radiance);
187     if (cloud_cover > 0) {
188     double zenithbr = get_zenith_brightness(sundir);
189     double skybr = get_overcast_brightness(sundir[2], zenithbr);
190     for (int i = 0; i < NSSAMP; ++i) {
191     sun_radiance[i] =
192     wmean2(sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
193     }
194     }
195    
196     FILE *rfp = fopen(radfile, "w");
197     write_rad_file(rfp, sun_radiance, sundir, skyfile, grndfile);
198     fclose(rfp);
199     return 1;
200     }
201    
202     static DpPaths get_dppaths(const double aod, const char *tag) {
203     DpPaths paths;
204    
205     snprintf(paths.tau, PATH_MAX, "tau_%s_%.2f.dat", tag, aod);
206     snprintf(paths.scat, PATH_MAX, "scat_%s_%.2f.dat", tag, aod);
207     snprintf(paths.scat1m, PATH_MAX, "scat1m_%s_%.2f.dat", tag, aod);
208     snprintf(paths.irrad, PATH_MAX, "irrad_%s_%.2f.dat", tag, aod);
209    
210     return paths;
211     }
212    
213     static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag, const int is_summer,
214     const double s_latitude) {
215     // Set rayleigh density profile
216     if (fabs(s_latitude*180.0 / PI) > ARCTIC_LAT) {
217     tag[0] = 's';
218     if (is_summer) {
219     tag[1] = 's';
220     atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SS;
221     atmos->beta_r0 = BR0_SS;
222     } else {
223     tag[1] = 'w';
224     atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_SW;
225     atmos->beta_r0 = BR0_SW;
226     }
227     } else if (fabs(s_latitude*180.0/PI) > TROPIC_LAT) {
228     tag[0] = 'm';
229     if (is_summer) {
230     tag[1] = 's';
231     atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MS;
232     atmos->beta_r0 = BR0_MS;
233     } else {
234     tag[1] = 'w';
235     atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_MW;
236     atmos->beta_r0 = BR0_MW;
237     }
238     } else {
239     tag[0] = 't';
240     tag[1] = 'r';
241     atmos->rayleigh_density.layers[0].exp_scale = -1.0 / HR_T;
242     atmos->beta_r0 = BR0_T;
243     }
244     tag[2] = '\0';
245     }
246    
247     static Atmosphere init_atmos(const double aod, const double grefl) {
248     Atmosphere atmos = {
249     .ozone_density = {.layers =
250     {
251     {.width = 25000.0,
252     .exp_term = 0.0,
253     .exp_scale = 0.0,
254     .linear_term = 1.0 / 15000.0,
255     .constant_term = -2.0 / 3.0},
256     {.width = AH,
257     .exp_term = 0.0,
258     .exp_scale = 0.0,
259     .linear_term = -1.0 / 15000.0,
260     .constant_term = 8.0 / 3.0},
261     }},
262     .rayleigh_density = {.layers =
263     {
264     {.width = AH,
265     .exp_term = 1.0,
266     .exp_scale = -1.0 / HR_MS,
267     .linear_term = 0.0,
268     .constant_term = 0.0},
269     }},
270     .beta_r0 = BR0_MS,
271     .beta_scale = aod / AOD0_CA,
272     .beta_m = NULL,
273     .grefl = grefl
274     };
275     return atmos;
276     }
277    
278     int main(int argc, char *argv[]) {
279     progname = argv[0];
280     int month, day;
281     double hour;
282     FVECT sundir;
283     int num_threads = 1;
284     int sorder = 4;
285     int year = 0;
286     int tsolar = 0;
287     double grefl = 0.2;
288     double ccover = 0.0;
289     int res = 128;
290     double aod = AOD0_CA;
291     char *outname = "out";
292     char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
293     char lstag[3];
294    
295     if (argc < 4) {
296     fprintf(stderr, "Usage: %s month day hour -y year -a lat -o lon -m tz -d aod -r res -n nproc -c ccover -l mie -g grefl -f outpath\n",
297     argv[0]);
298     return 0;
299     }
300    
301     month = atoi(argv[1]);
302     day = atoi(argv[2]);
303     hour = atof(argv[3]);
304    
305     if (!compute_sundir(year, month, day, hour, tsolar, sundir)) {
306     fprintf(stderr, "Cannot compute solar angle\n");
307     exit(1);
308     }
309    
310     for (int i = 4; i < argc; i++) {
311     if (argv[i][0] == '-') {
312     switch (argv[i][1]) {
313     case 'a':
314     s_latitude = atof(argv[++i]) * (PI / 180.0);
315     break;
316     case 'g':
317     grefl = atof(argv[++i]);
318     break;
319     case 'c':
320     ccover = atof(argv[++i]);
321     break;
322     case 'd':
323     aod = atof(argv[++i]);
324     break;
325     case 'i':
326     sorder = atoi(argv[++i]);
327     break;
328     case 'l':
329     mie_path = argv[++i];
330     break;
331     case 'm':
332     s_meridian = atof(argv[++i]) * (PI / 180.0);
333     break;
334     case 'o':
335     s_longitude = atof(argv[++i]) * (PI / 180.0);
336     break;
337     case 'n':
338     num_threads = atoi(argv[++i]);
339     break;
340     case 'y':
341     year = atoi(argv[++i]);
342     break;
343     case 'f':
344     outname = argv[++i];
345     break;
346     case 'r':
347     res = atoi(argv[++i]);
348     break;
349     default:
350     fprintf(stderr, "Unknown option %s\n", argv[i]);
351     exit(1);
352     }
353     }
354     }
355    
356     Atmosphere clear_atmos = init_atmos(aod, grefl);
357    
358     int is_summer = (month >= SUMMER_START && month <= SUMMER_END);
359     if (s_latitude < 0) {
360     is_summer = !is_summer;
361     }
362     set_rayleigh_density_profile(&clear_atmos, lstag, is_summer, s_latitude);
363    
364     // Load mie density data
365     DATARRAY *mie_dp = getdata(mie_path);
366     if (mie_dp == NULL) {
367     fprintf(stderr, "Error reading mie data\n");
368     return 0;
369     }
370     clear_atmos.beta_m = mie_dp;
371    
372     DpPaths clear_paths = get_dppaths(aod, lstag);
373    
374     if (getpath(clear_paths.tau, ".", R_OK) == NULL ||
375     getpath(clear_paths.scat, ".", R_OK) == NULL ||
376     getpath(clear_paths.scat1m, ".", R_OK) == NULL ||
377     getpath(clear_paths.irrad, ".", R_OK) == NULL) {
378     printf("# Precomputing...\n");
379     if (!precompute(sorder, clear_paths, &clear_atmos, num_threads)) {
380     fprintf(stderr, "Precompute failed\n");
381     return 0;
382     }
383     }
384    
385     DATARRAY *tau_clear_dp = getdata(clear_paths.tau);
386     DATARRAY *irrad_clear_dp = getdata(clear_paths.irrad);
387     DATARRAY *scat_clear_dp = getdata(clear_paths.scat);
388     DATARRAY *scat1m_clear_dp = getdata(clear_paths.scat1m);
389    
390     if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
391     irrad_clear_dp, ccover, sundir, grefl, res, outname)) {
392     fprintf(stderr, "gen_spect_sky failed\n");
393     exit(1);
394     }
395    
396     freedata(mie_dp);
397     freedata(tau_clear_dp);
398     freedata(scat_clear_dp);
399     freedata(irrad_clear_dp);
400     freedata(scat1m_clear_dp);
401    
402     return 1;
403     }