ViewVC Help

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

root/Development/ray/src/gen/genssky.c

Comparing ray/src/gen/genssky.c (file contents):
Revision 2.2 by greg, Fri Jul 19 23:38:28 2024 UTC vs.
Revision 2.4 by greg, Thu Aug 8 02:00:48 2024 UTC

#	Line 8 | Line 8 | static const char RCSid[] = "$Id$";
8		#include "copyright.h"
9		#include "resolu.h"
10		#include "rtio.h"
11	–	#include "view.h"
11		#include <ctype.h>
12		#ifdef _WIN32
13		#include <windows.h>
#	Line 61 | Line 60 | static int make_directory(const char *path) {
60		#endif
61		}
62
63	+	inline static float deg2rad(float deg) { return deg * (PI / 180.); }
64	+
65		static int cvthour(char *hs, int *tsolar, double *hour) {
66		char *cp = hs;
67		int i, j;
#	Line 131 | Line 132 | static void basename(const char *path, char *output, s
132		}
133		}
134
135	<	char *join_paths(const char *path1, const char *path2) {
135	>	static char *join_paths(const char *path1, const char *path2) {
136		size_t len1 = strlen(path1);
137		size_t len2 = strlen(path2);
138		int need_separator = (path1[len1 - 1] != DIRSEP);
#	Line 178 | Line 179 | static double get_overcast_brightness(const double dz,
179
180		static void write_header(const int argc, char **argv, const double cloud_cover,
181		const double grefl, const int res) {
182	+	int i;
183		printf("# ");
184	<	for (int i = 0; i < argc; i++) {
184	>	for (i = 0; i < argc; i++) {
185		printf("%s ", argv[i]);
186		}
187		printf("\n");
188	<	printf("#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
189	<	cloud_cover, grefl, res);
188	>	printf(
189	>	"#Cloud cover: %g\n#Ground reflectance: %g\n#Sky map resolution: %d\n\n",
190	>	cloud_cover, grefl, res);
191		}
192
193		static void write_rad(const double *sun_radiance, const FVECT sundir,
194	<	const char skyfile[PATH_MAX],
192	<	const char grndfile[PATH_MAX]) {
194	>	const char *ddir, const char *skyfile) {
195		if (sundir[2] > 0) {
196		printf("void spectrum sunrad\n0\n0\n22 380 780 ");
197		/* Normalize to one */
198		double sum = 0.0;
199	<	for (int i = 0; i < NSSAMP; ++i) {
199	>	int i;
200	>	for (i = 0; i < NSSAMP; ++i) {
201		sum += sun_radiance[i];
202		}
203		double mean = sum / NSSAMP;
204	<	for (int i = 0; i < NSSAMP; ++i) {
204	>	for (i = 0; i < NSSAMP; ++i) {
205		printf("%.3f ", sun_radiance[i] / mean);
206		}
207		double intensity = mean * WVLSPAN;
#	Line 207 | Line 210 | static void write_rad(const double *sun_radiance, cons
210		printf("solar source sun\n0\n0\n4 %f %f %f 0.533\n\n", sundir[0], sundir[1],
211		sundir[2]);
212		}
213	<	printf("void specpict skymap\n8 noop %s fisheye.cal fish_u fish_v -rx 90 "
214	<	"-mx\n0\n0\n\n",
213	>	printf("void specpict skyfunc\n5 noop %s . 'atan2(Dy,Dx)/PI+1' "
214	>	"'acos(Dz)/PI'\n0\n0\n\n",
215		skyfile);
213	–
214	–	printf("void specpict grndmap\n8 noop %s fisheye.cal fish_u fish_v -rx -90 "
215	–	"-my\n0\n0\n\n",
216	–	grndfile);
217	–	printf("void mixfunc skyfunc\n4 skymap grndmap if(Dz,1,0) .\n0\n0\n");
216		}
217
218		static void write_hsr_header(FILE *fp, RESOLU *res) {
219	<	float wvsplit[4] = {380, 480, 588,
222	<	780}; /* RGB wavelength limits+partitions (nm) */
219	>	float wvsplit[4] = {380, 480, 588, 780};
220		newheader("RADIANCE", fp);
221		fputncomp(NSSAMP, fp);
222		fputwlsplit(wvsplit, fp);
#	Line 228 | Line 225 | static void write_hsr_header(FILE *fp, RESOLU *res) {
225		fputsresolu(res, fp);
226		}
227
228	+	static inline float frac(float x) { return x - floor(x); }
229	+
230		int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_clear,
231		DATARRAY *scat1m_clear, DATARRAY *irrad_clear,
232		const double cloud_cover, const FVECT sundir,
233	<	const double grefl, const int res, const char *outname) {
233	>	const double grefl, const int res, const char *outname,
234	>	const char *ddir) {
235		char skyfile[PATH_MAX];
236		char grndfile[PATH_MAX];
237	<	if (!snprintf(skyfile, sizeof(skyfile), "%s_sky.hsr", outname)) {
237	>	if (!snprintf(skyfile, sizeof(skyfile), "%s%c%s_sky.hsr", ddir, DIRSEP,
238	>	outname)) {
239		fprintf(stderr, "Error setting sky file name\n");
240		return 0;
241		};
242	<	if (!snprintf(grndfile, sizeof(grndfile), "%s_ground.hsr", outname)) {
243	<	fprintf(stderr, "Error setting ground file name\n");
244	<	return 0;
244	<	}
245	<	RESOLU rs = {PIXSTANDARD, res, res};
242	>	int xres = res;
243	>	int yres = xres / 2;
244	>	RESOLU rs = {PIXSTANDARD, xres, yres};
245		FILE *skyfp = fopen(skyfile, "w");
247	–	FILE *grndfp = fopen(grndfile, "w");
248	–	write_hsr_header(grndfp, &rs);
246		write_hsr_header(skyfp, &rs);
250	–	VIEW skyview = {VT_ANG, {0., 0., 0.}, {0., 0., 1.}, {0., 1., 0.}, 1.,
251	–	180.,   180.,         0.,           0.,           0.,
252	–	0.,     {0., 0., 0.}, {0., 0., 0.}, 0.,           0.};
253	–	VIEW grndview = {
254	–	VT_ANG, {0., 0., 0.}, {0., 0., -1.}, {0., 1., 0.}, 1., 180., 180., 0., 0.,
255	–	0.,     0.,           {0., 0., 0.},  {0., 0., 0.}, 0., 0.};
256	–	setview(&skyview);
257	–	setview(&grndview);
247
248		CNDX[3] = NSSAMP;
249
250	<	FVECT view_point = {0, 0, ER};
250	>	FVECT view_point = {0, 0, ER + 10};
251		const double radius = VLEN(view_point);
252		const double sun_ct = fdot(view_point, sundir) / radius;
253	<	for (unsigned int j = 0; j < res; ++j) {
254	<	for (unsigned int i = 0; i < res; ++i) {
255	<	RREAL loc[2];
256	<	FVECT rorg = {0};
268	<	FVECT rdir_sky = {0};
269	<	FVECT rdir_grnd = {0};
270	<	SCOLOR sky_radiance = {0};
271	<	SCOLOR ground_radiance = {0};
253	>	int i, j, k;
254	>	for (j = 0; j < yres; ++j) {
255	>	for (i = 0; i < xres; ++i) {
256	>	SCOLOR radiance = {0};
257		SCOLR sky_sclr = {0};
273	–	SCOLR ground_sclr = {0};
258
259	<	pix2loc(loc, &rs, i, j);
260	<	viewray(rorg, rdir_sky, &skyview, loc[0], loc[1]);
261	<	viewray(rorg, rdir_grnd, &grndview, loc[0], loc[1]);
259	>	float px = i / (xres - 1.0);
260	>	float py = j / (yres - 1.0);
261	>	float lambda = ((1 - py) * PI) - (PI / 2.0);
262	>	float phi = (px * 2.0 * PI) - PI;
263
264	<	const double mu_sky = fdot(view_point, rdir_sky) / radius;
265	<	const double nu_sky = fdot(rdir_sky, sundir);
264	>	FVECT rdir = {cos(lambda) * cos(phi), cos(lambda) * sin(phi),
265	>	sin(lambda)};
266
267	<	const double mu_grnd = fdot(view_point, rdir_grnd) / radius;
268	<	const double nu_grnd = fdot(rdir_grnd, sundir);
267	>	const double mu = fdot(view_point, rdir) / radius;
268	>	const double nu = fdot(rdir, sundir);
269
270	<	get_sky_radiance(scat_clear, scat1m_clear, radius, mu_sky, sun_ct, nu_sky,
271	<	sky_radiance);
272	<	get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
273	<	view_point, rdir_grnd, radius, mu_grnd, sun_ct,
274	<	nu_grnd, grefl, sundir, ground_radiance);
270	>	/* hit ground */
271	>	if (rdir[2] < 0) {
272	>	get_ground_radiance(tau_clear, scat_clear, scat1m_clear, irrad_clear,
273	>	view_point, rdir, radius, mu, sun_ct, nu, grefl,
274	>	sundir, radiance);
275	>	} else {
276	>	get_sky_radiance(scat_clear, scat1m_clear, radius, mu, sun_ct, nu,
277	>	radiance);
278	>	}
279
280	<	for (int k = 0; k < NSSAMP; ++k) {
281	<	sky_radiance[k] *= WVLSPAN;
293	<	ground_radiance[k] *= WVLSPAN;
280	>	for (k = 0; k < NSSAMP; ++k) {
281	>	radiance[k] *= WVLSPAN;
282		}
283
284		if (cloud_cover > 0) {
285		double zenithbr = get_zenith_brightness(sundir);
286		double grndbr = zenithbr * GNORM;
287	<	double skybr = get_overcast_brightness(rdir_sky[2], zenithbr);
288	<	for (int k = 0; k < NSSAMP; ++k) {
289	<	sky_radiance[k] =
290	<	wmean2(sky_radiance[k], skybr * D6415[k], cloud_cover);
291	<	ground_radiance[k] =
292	<	wmean2(ground_radiance[k], grndbr * D6415[k], cloud_cover);
287	>	double skybr = get_overcast_brightness(rdir[2], zenithbr);
288	>	if (rdir[2] < 0) {
289	>	for (k = 0; k < NSSAMP; ++k) {
290	>	radiance[k] = wmean2(radiance[k], grndbr * D6415[k], cloud_cover);
291	>	}
292	>	} else {
293	>	for (k = 0; k < NSSAMP; ++k) {
294	>	radiance[k] = wmean2(radiance[k], skybr * D6415[k], cloud_cover);
295	>	}
296		}
297		}
298
299	<	scolor2scolr(sky_sclr, sky_radiance, 20);
299	>	scolor2scolr(sky_sclr, radiance, 20);
300		putbinary(sky_sclr, LSCOLR, 1, skyfp);
310	–
311	–	scolor2scolr(ground_sclr, ground_radiance, 20);
312	–	putbinary(ground_sclr, LSCOLR, 1, grndfp);
301		}
302		}
303		fclose(skyfp);
316	–	fclose(grndfp);
304
305		/* Get solar radiance */
306		double sun_radiance[NSSAMP] = {0};
#	Line 322 | Line 309 | int gen_spect_sky(DATARRAY *tau_clear, DATARRAY *scat_
309		if (cloud_cover > 0) {
310		double zenithbr = get_zenith_brightness(sundir);
311		double skybr = get_overcast_brightness(sundir[2], zenithbr);
312	<	for (int i = 0; i < NSSAMP; ++i) {
312	>	int i;
313	>	for (i = 0; i < NSSAMP; ++i) {
314		sun_radiance[i] =
315		wmean2(sun_radiance[i], D6415[i] * skybr / WVLSPAN, cloud_cover);
316		}
317		}
318
319	<	write_rad(sun_radiance, sundir, skyfile, grndfile);
319	>	write_rad(sun_radiance, sundir, ddir, skyfile);
320		return 1;
321		}
322
#	Line 351 | Line 339 | static DpPaths get_dppaths(const char *dir, const doub
339		static void set_rayleigh_density_profile(Atmosphere *atmos, char *tag,
340		const int is_summer,
341		const double s_latitude) {
354	–	/* Set rayleigh density profile */
342		if (fabs(s_latitude * 180.0 / PI) > ARCTIC_LAT) {
343		tag[0] = 's';
344		if (is_summer) {
#	Line 424 | Line 411 | int main(int argc, char *argv[]) {
411		int got_meridian = 0;
412		double grefl = 0.2;
413		double ccover = 0.0;
414	<	int res = 128;
414	>	int res = 64;
415		double aod = AOD0_CA;
416		char *outname = "out";
417		char *mie_path = getpath("mie_ca.dat", getrlibpath(), R_OK);
418		char mie_name[20] = "mie_ca";
419		char lstag[3];
420		char *ddir = ".";
421	+	int i;
422
423	<	if (!strcmp(argv[1], "-defaults")) {
423	>	if (argc == 2 && !strcmp(argv[1], "-defaults")) {
424		printf("-i %d\t\t\t\t#scattering order\n", sorder);
425		printf("-g %f\t\t\t#ground reflectance\n", grefl);
426		printf("-c %f\t\t\t#cloud cover\n", ccover);
#	Line 440 | Line 428 | int main(int argc, char *argv[]) {
428		printf("-d %f\t\t\t#broadband aerosol optical depth\n", AOD0_CA);
429		printf("-f %s\t\t\t\t#output name (-f)\n", outname);
430		printf("-p %s\t\t\t\t#atmos data directory\n", ddir);
431	<	exit(1);
431	>	exit(0);
432		}
433
434		if (argc < 4) {
#	Line 468 | Line 456 | int main(int argc, char *argv[]) {
456		exit(1);
457		}
458
459	<	for (int i = 4; i < argc; i++) {
459	>	for (i = 4; i < argc; i++) {
460		if (argv[i][0] == '-') {
461		switch (argv[i][1]) {
462		case 'a':
#	Line 547 | Line 535 | int main(int argc, char *argv[]) {
535
536		char gsdir[PATH_MAX];
537		size_t siz = strlen(ddir);
538	<	if (ISDIRSEP(ddir[siz-1]))
539	<	ddir[siz-1] = '\0';
538	>	if (ISDIRSEP(ddir[siz - 1]))
539	>	ddir[siz - 1] = '\0';
540		snprintf(gsdir, PATH_MAX, "%s%catmos_data", ddir, DIRSEP);
553	–	printf("gsdir: %s\n", gsdir);
541		if (!make_directory(gsdir)) {
542		fprintf(stderr, "Failed creating atmos_data directory");
543		exit(1);
#	Line 576 | Line 563 | int main(int argc, char *argv[]) {
563		write_header(argc, argv, ccover, grefl, res);
564
565		if (!gen_spect_sky(tau_clear_dp, scat_clear_dp, scat1m_clear_dp,
566	<	irrad_clear_dp, ccover, sundir, grefl, res, outname)) {
566	>	irrad_clear_dp, ccover, sundir, grefl, res, outname,
567	>	ddir)) {
568		fprintf(stderr, "gen_spect_sky failed\n");
569		exit(1);
570		}

Diff Legend

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