Output: write HDR FITS with camera metadata

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)
 

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

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;

src/optics.c

@@ -8,6 +8,10 @@
 #include <stdlib.h>
 #include <string.h>
 
+#ifdef ENABLE_HDR_DEBUG
+#include <fitsio.h>
+#endif
+
 #ifdef ENABLE_PNG
 #include <png.h>
 #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

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