From d5d21c7659fb9019c339bd9947e875aa24294ab4 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Fri, 28 Aug 2026 19:56:07 -0400 Subject: [PATCH] Output: write HDR FITS with camera metadata --- Makefile | 4 +-- README.md | 11 ++++-- src/main.c | 5 +-- src/optics.c | 95 +++++++++++++++++++++++++++++++++++++++++----------- src/optics.h | 5 +-- 5 files changed, 92 insertions(+), 28 deletions(-) diff --git a/Makefile b/Makefile index 5f52ada..2f8bdb5 100644 --- a/Makefile +++ b/Makefile @@ -52,9 +52,9 @@ minkowski: schwarzschild: $(MAKE) SPACETIME=schwarzschild all -# Deliberately separate from the production binary: enables --hdr-output PFM. +# Deliberately separate from the production binary: enables --hdr-output FITS. $(PSF_HDR_TEST_TARGET): $(COMMON_SOURCES) src/spacetime_minkowski.c src/main.c | build - $(CC) $(CPPFLAGS) -DENABLE_HDR_DEBUG -DSPACETIME_MINKOWSKI $(CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ + $(CC) $(CPPFLAGS) -DENABLE_HDR_DEBUG -DSPACETIME_MINKOWSKI $(CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) $(shell pkg-config --libs cfitsio) -o $@ psf-hdr-test: $(PSF_HDR_TEST_TARGET) diff --git a/README.md b/README.md index 5d41a63..e67cf57 100644 --- a/README.md +++ b/README.md @@ -174,8 +174,15 @@ the ray counts before each slab is loaded. For PSF validation only, `make psf-hdr-test` builds `build/minkowski_psf_hdr_test`, a separate binary with a `--hdr-output PATH` -option. It writes the pre-tone-mapping RGB framebuffer as a 32-bit float PFM -image; the ordinary binaries do not contain this option or writer. +option. It writes the pre-tone-mapping RGB framebuffer as a three-plane, +32-bit float FITS image. Values remain linear HDR at the renderer's arbitrary +scale; no tone mapping or per-frame normalization is applied. The ordinary +binaries do not contain this option or writer. + +The FITS header describes a synthetic 8640-by-5760, 36-by-24 mm full-frame +sensor with 4.1667 um pixels. Each render records its active centered crop and +derives `FOCALLEN` from that render's width and horizontal `--fov-deg`; these +camera fields support plate-solving workflows but do not calibrate flux. Pass `--draw-mesh` to alpha-composite image-plane triangle edges as one-pixel-wide 0.5 linear-gray diagnostic lines at 0.5 opacity. The line diff --git a/src/main.c b/src/main.c index 1df5a91..5551e7d 100644 --- a/src/main.c +++ b/src/main.c @@ -312,8 +312,9 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, frame_draw_mesh(&mesh, hdr, s->width, s->height, 0.5, 0.5); #ifdef ENABLE_HDR_DEBUG if (s->hdr_output_path != NULL && - write_hdr_pfm(s->hdr_output_path, hdr, s->width, s->height)) { - perror(s->hdr_output_path); + write_hdr_fits(s->hdr_output_path, hdr, s->width, s->height, + s->horizontal_fov_deg)) { + fprintf(stderr, "Failed to write HDR FITS image: %s\n", s->hdr_output_path); frame_lens_mesh_destroy(&mesh); free(hdr); return -1; diff --git a/src/optics.c b/src/optics.c index 5ac7fcb..32dccb8 100644 --- a/src/optics.c +++ b/src/optics.c @@ -8,6 +8,10 @@ #include #include +#ifdef ENABLE_HDR_DEBUG +#include +#endif + #ifdef ENABLE_PNG #include #endif @@ -446,28 +450,79 @@ int write_tonemapped_image(const char *path, const double *hdr, int width, int h } #ifdef ENABLE_HDR_DEBUG -int write_hdr_pfm(const char *path, const double *hdr, int width, int height) +int write_hdr_fits(const char *path, const double *hdr, int width, int height, + double horizontal_fov_deg) { - if (path == NULL || hdr == NULL || width <= 0 || height <= 0) + if (path == NULL || hdr == NULL || width <= 0 || height <= 0 || + !isfinite(horizontal_fov_deg) || horizontal_fov_deg <= 0.0 || + horizontal_fov_deg >= 179.0) return -1; - FILE *file = fopen(path, "wb"); - if (file == NULL) - return -1; - const uint16_t endian_probe = 1; - const char *scale = *(const unsigned char *)&endian_probe == 1 ? "-1.0" : "1.0"; - int result = fprintf(file, "PF\n%d %d\n%s\n", width, height, scale) < 0 ? -1 : 0; - /* PFM rows are stored bottom-to-top. Its negative scale declares little-endian - * float samples, avoiding an unnecessary byte swap on the normal test host. */ - for (int row = height - 1; result == 0 && row >= 0; --row) - for (int column = 0; column < width * 3; ++column) { - const float sample = (float)hdr[(size_t)row * width * 3 + column]; - if (fwrite(&sample, sizeof sample, 1, file) != 1) { - result = -1; - break; - } + const size_t path_length = strlen(path); + char *overwrite_path = malloc(path_length + 2); + float *scanline = malloc((size_t)width * sizeof *scanline); + fitsfile *file = NULL; + int status = 0; + int result = -1; + long dimensions[3] = {width, height, 3}; + char creator[] = "GR_4d_raytracing"; + char unit[] = "linear HDR; arbitrary renderer scale"; + char color_axis[] = "RGB"; + char instrument[] = "GR4D virtual full-frame 50MP"; + char detector_size[] = "[1:8640,1:5760]"; + char crop_size[FLEN_VALUE]; + double pixel_size_um = 36e3 / 8640.0; + const double crop_width_mm = width * pixel_size_um / 1000.0; + const double crop_height_mm = height * pixel_size_um / 1000.0; + double focal_length_mm = + crop_width_mm / (2.0 * tan(horizontal_fov_deg * pi / 360.0)); + if (overwrite_path == NULL || scanline == NULL) + goto done; + snprintf(crop_size, sizeof crop_size, "%.6g x %.6g mm centered crop", + crop_width_mm, crop_height_mm); + /* CFITSIO uses a leading ! to give this output the same overwrite behavior + * as the former PFM writer. */ + overwrite_path[0] = '!'; + memcpy(overwrite_path + 1, path, path_length + 1); + fits_create_file(&file, overwrite_path, &status); + fits_create_img(file, FLOAT_IMG, 3, dimensions, &status); + fits_update_key(file, TSTRING, "CREATOR", creator, NULL, &status); + fits_update_key(file, TSTRING, "BUNIT", unit, NULL, &status); + fits_update_key(file, TSTRING, "CTYPE3", color_axis, NULL, &status); + fits_update_key(file, TSTRING, "INSTRUME", instrument, + "synthetic camera metadata", &status); + fits_update_key(file, TSTRING, "DETSIZE", detector_size, + "full-frame detector pixels", &status); + fits_update_key(file, TDOUBLE, "FOCALLEN", &focal_length_mm, + "mm; derived from horizontal field of view", &status); + fits_update_key(file, TDOUBLE, "XPIXSZ", &pixel_size_um, + "um; full-frame detector pixel size", &status); + fits_update_key(file, TDOUBLE, "YPIXSZ", &pixel_size_um, + "um; full-frame detector pixel size", &status); + fits_update_key(file, TSTRING, "CROPSIZE", crop_size, + "active sensor area for this render", &status); + /* FITS image coordinates start at the lower-left. The HDR buffer is + * top-down like PNG, so reverse rows while preserving its visual orientation. */ + for (int channel = 0; status == 0 && channel < 3; ++channel) + for (int row = 0; status == 0 && row < height; ++row) { + const size_t offset = (size_t)row * width * 3; + for (int column = 0; column < width; ++column) + scanline[column] = (float)hdr[offset + 3 * column + channel]; + long first_pixel[3] = {1, height - row, channel + 1}; + fits_write_pix(file, TFLOAT, first_pixel, width, scanline, &status); } - if (fclose(file) != 0) - result = -1; - return result; + if (status == 0) + result = 0; +done: + if (file != NULL) { + int close_status = 0; + fits_close_file(file, &close_status); + if (status == 0 && close_status != 0) + status = close_status; + } + if (status != 0) + fits_report_error(stderr, status); + free(scanline); + free(overwrite_path); + return status == 0 ? result : -1; } #endif diff --git a/src/optics.h b/src/optics.h index 08a53f9..c8ee4f8 100644 --- a/src/optics.h +++ b/src/optics.h @@ -51,8 +51,9 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y, int write_tonemapped_image(const char *path, const double *hdr, int width, int height); #ifdef ENABLE_HDR_DEBUG -/* Test-build-only: writes the pre-tone-mapping framebuffer as RGB float PFM. */ -int write_hdr_pfm(const char *path, const double *hdr, int width, int height); +/* Test-build-only: writes the pre-tone-mapping framebuffer as 32-bit RGB FITS. */ +int write_hdr_fits(const char *path, const double *hdr, int width, int height, + double horizontal_fov_deg); #endif #endif