Add lazy tiled 2MASS catalog support

src/catalog.c

@@ -1,6 +1,7 @@
 #include "catalog.h"
 
 #include <math.h>
+#include <limits.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -65,8 +66,11 @@ int catalog_load_csv(StarCatalog *catalog, const char *path)
     FILE *file = fopen(path, "r");
     char line[256];
     size_t capacity = 0;
-    catalog->stars = NULL;
-    catalog->count = 0;
+    if (catalog == NULL) {
+        if (file != NULL) fclose(file);
+        return -1;
+    }
+    *catalog = (StarCatalog){0};
     if (file == NULL || fgets(line, sizeof line, file) == NULL) goto fail;
 
     while (fgets(line, sizeof line, file) != NULL) {
@@ -99,9 +103,238 @@ fail:
     return -1;
 }
 
+static size_t tile_index(int ra_index, int dec_index)
+{
+    return (size_t)dec_index * CATALOG_ALL_SKY_RA_TILES + ra_index;
+}
+
+static int load_tile(StarCatalog *catalog, int ra_index, int dec_index)
+{
+    CatalogTile *tile = &catalog->tiles[tile_index(ra_index, dec_index)];
+    char path[PATH_MAX];
+    int written;
+    if (tile->state != 0)
+        return tile->state == 1 ? 0 : -1;
+    written = snprintf(path, sizeof path, "%s/tile_ra%03d_dec%03d.csv",
+                       catalog->all_sky_root, ra_index, dec_index);
+    if (written < 0 || (size_t)written >= sizeof path) {
+        tile->state = -1;
+        return -1;
+    }
+    StarCatalog temporary = {0};
+    if (catalog_load_csv(&temporary, path)) {
+        tile->state = -1;
+        return -1;
+    }
+    tile->stars = temporary.stars;
+    tile->count = temporary.count;
+    tile->state = 1;
+    catalog->count += tile->count;
+    return 0;
+}
+
+int catalog_load_all_sky(StarCatalog *catalog, const char *directory)
+{
+    size_t tile_count = (size_t)CATALOG_ALL_SKY_RA_TILES *
+                        CATALOG_ALL_SKY_DEC_TILES;
+    if (catalog == NULL || directory == NULL || directory[0] == '\0')
+        return -1;
+    *catalog = (StarCatalog){0};
+    catalog->all_sky_root = malloc(strlen(directory) + 1);
+    catalog->tiles = calloc(tile_count, sizeof *catalog->tiles);
+    if (catalog->all_sky_root == NULL || catalog->tiles == NULL) {
+        catalog_destroy(catalog);
+        return -1;
+    }
+    strcpy(catalog->all_sky_root, directory);
+    catalog->kind = STAR_CATALOG_ALL_SKY;
+    return 0;
+}
+
+static void lon_lat_from_direction(const double direction[3], double *longitude,
+                                   double *latitude)
+{
+    *longitude = atan2(direction[2], direction[0]) * 180.0 / PI;
+    if (*longitude < 0.0)
+        *longitude += 360.0;
+    *latitude = asin(fmax(-1.0, fmin(1.0, direction[1]))) * 180.0 / PI;
+}
+
+static double dot(const double a[3], const double b[3])
+{
+    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+
+static void cross(const double a[3], const double b[3], double out[3])
+{
+    out[0] = a[1] * b[2] - a[2] * b[1];
+    out[1] = a[2] * b[0] - a[0] * b[2];
+    out[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+static void direction_from_lon_lat(double longitude, double latitude,
+                                   double direction[3])
+{
+    const double lon = longitude * PI / 180.0;
+    const double lat = latitude * PI / 180.0;
+    const double cos_lat = cos(lat);
+    direction[0] = cos_lat * cos(lon);
+    direction[1] = sin(lat);
+    direction[2] = cos_lat * sin(lon);
+}
+
+/* The minimum of a plane dot-product over a longitude/latitude rectangle is
+ * attained on a boundary or at its antipodal stationary point.  This gives a
+ * conservative, analytic whole-tile containment test; it is not a corner-only
+ * approximation. */
+static double tile_plane_minimum(const double normal[3], double lon_lo,
+                                 double lon_hi, double lat_lo, double lat_hi)
+{
+    double values[32];
+    size_t count = 0;
+    const double phase = atan2(normal[2], normal[0]) * 180.0 / PI;
+    const double latitude_phase = atan2(normal[1], hypot(normal[0], normal[2])) * 180.0 / PI;
+    const double candidate_lon[] = {lon_lo, lon_hi, phase + 180.0, phase - 180.0};
+    const double candidate_lat[] = {lat_lo, lat_hi, latitude_phase + 180.0,
+                                    latitude_phase - 180.0};
+    for (size_t i = 0; i < sizeof candidate_lon / sizeof *candidate_lon; ++i)
+        for (size_t j = 0; j < 2; ++j) {
+            double lon = candidate_lon[i];
+            while (lon < lon_lo) lon += 360.0;
+            while (lon > lon_hi) lon -= 360.0;
+            if (lon >= lon_lo && lon <= lon_hi) {
+                double point[3];
+                direction_from_lon_lat(lon, j == 0 ? lat_lo : lat_hi, point);
+                values[count++] = dot(normal, point);
+            }
+        }
+    for (size_t i = 0; i < sizeof candidate_lon / sizeof *candidate_lon; ++i)
+        for (size_t j = 0; j < sizeof candidate_lat / sizeof *candidate_lat; ++j) {
+            double lon = candidate_lon[i];
+            while (lon < lon_lo) lon += 360.0;
+            while (lon > lon_hi) lon -= 360.0;
+            if (lon >= lon_lo && lon <= lon_hi &&
+                candidate_lat[j] >= lat_lo && candidate_lat[j] <= lat_hi) {
+                double point[3];
+                direction_from_lon_lat(lon, candidate_lat[j], point);
+                values[count++] = dot(normal, point);
+            }
+        }
+    for (size_t i = 0; i < 2; ++i)
+        for (size_t j = 0; j < sizeof candidate_lat / sizeof *candidate_lat; ++j)
+            if (candidate_lat[j] >= lat_lo && candidate_lat[j] <= lat_hi) {
+                double point[3];
+                direction_from_lon_lat(i == 0 ? lon_lo : lon_hi,
+                                       candidate_lat[j], point);
+                values[count++] = dot(normal, point);
+            }
+    double minimum = values[0];
+    for (size_t i = 1; i < count; ++i)
+        if (values[i] < minimum) minimum = values[i];
+    return minimum;
+}
+
+static int tile_is_fully_contained(const double direction[3][3], int ra_index,
+                                   int dec_index)
+{
+    const double lon_lo = ra_index;
+    const double lon_hi = ra_index + 1.0;
+    const double lat_lo = dec_index - 90.0;
+    const double lat_hi = lat_lo + 1.0;
+    for (int edge = 0; edge < 3; ++edge) {
+        const double *left = direction[edge];
+        const double *right = direction[(edge + 1) % 3];
+        const double *opposite = direction[(edge + 2) % 3];
+        double normal[3];
+        cross(left, right, normal);
+        if (dot(normal, opposite) < 0.0)
+            for (int axis = 0; axis < 3; ++axis) normal[axis] = -normal[axis];
+        if (tile_plane_minimum(normal, lon_lo, lon_hi, lat_lo, lat_hi) < -1e-14)
+            return 0;
+    }
+    return 1;
+}
+
+int catalog_visit_source_triangle(StarCatalog *catalog,
+                                  const double direction[3][3],
+                                  int load_missing, CatalogTileVisitor visitor,
+                                  void *context)
+{
+    if (catalog == NULL || direction == NULL || visitor == NULL)
+        return -1;
+    if (catalog->kind == STAR_CATALOG_MEMORY)
+        return visitor(catalog->stars, catalog->count, 0, context);
+    double longitude[3], latitude[3], unwrapped[3];
+    for (int i = 0; i < 3; ++i) {
+        lon_lat_from_direction(direction[i], &longitude[i], &latitude[i]);
+        unwrapped[i] = longitude[i];
+        while (unwrapped[i] - longitude[0] > 180.0) unwrapped[i] -= 360.0;
+        while (unwrapped[i] - longitude[0] < -180.0) unwrapped[i] += 360.0;
+    }
+    double lon_min = unwrapped[0], lon_max = unwrapped[0];
+    double lat_min = latitude[0], lat_max = latitude[0];
+    for (int i = 1; i < 3; ++i) {
+        if (unwrapped[i] < lon_min) lon_min = unwrapped[i];
+        if (unwrapped[i] > lon_max) lon_max = unwrapped[i];
+        if (latitude[i] < lat_min) lat_min = latitude[i];
+        if (latitude[i] > lat_max) lat_max = latitude[i];
+    }
+    /* A triangle containing a pole covers every RA there. */
+    const double north[3] = {0.0, 1.0, 0.0};
+    const double south[3] = {0.0, -1.0, 0.0};
+    int all_ra = 0;
+    for (int pole = 0; pole < 2; ++pole) {
+        const double *point = pole == 0 ? north : south;
+        int inside = 1;
+        for (int edge = 0; edge < 3; ++edge) {
+            double normal[3];
+            cross(direction[edge], direction[(edge + 1) % 3], normal);
+            if (dot(normal, point) * dot(normal, direction[(edge + 2) % 3]) < -1e-14) inside = 0;
+        }
+        if (inside) all_ra = 1;
+    }
+    /* Source edges are great-circle arcs, so their RA/Dec extrema need not be
+     * vertices.  This small guard band covers that curvature without pulling
+     * in an otherwise unrelated one-degree tile ring. */
+    lon_min -= 0.01;
+    lon_max += 0.01;
+    lat_min -= 0.01;
+    lat_max += 0.01;
+    if (lon_max - lon_min >= 360.0)
+        all_ra = 1;
+    const int dec_first = fmax(0, (int)floor(lat_min + 90.0));
+    const int dec_last = fmin(CATALOG_ALL_SKY_DEC_TILES - 1,
+                              (int)floor(lat_max + 90.0));
+    const int ra_first = (int)floor(lon_min);
+    const int ra_last = (int)floor(lon_max);
+    for (int dec = dec_first; dec <= dec_last; ++dec)
+        for (int raw_ra = all_ra ? 0 : ra_first;
+             raw_ra <= (all_ra ? 359 : ra_last); ++raw_ra) {
+            const int ra = (raw_ra % 360 + 360) % 360;
+            CatalogTile *tile = &catalog->tiles[tile_index(ra, dec)];
+            if (load_missing && load_tile(catalog, ra, dec))
+                continue; /* Downloader has not finished this tile yet. */
+            if (tile->state != 1 ||
+                visitor(tile->stars, tile->count,
+                        tile_is_fully_contained(direction, ra, dec), context))
+                return -1;
+        }
+    return 0;
+}
+
 void catalog_destroy(StarCatalog *catalog)
 {
+    if (catalog->tiles != NULL)
+        for (size_t i = 0;
+             i < (size_t)CATALOG_ALL_SKY_RA_TILES * CATALOG_ALL_SKY_DEC_TILES;
+             ++i)
+            free(catalog->tiles[i].stars);
     free(catalog->stars);
+    free(catalog->tiles);
+    free(catalog->all_sky_root);
     catalog->stars = NULL;
     catalog->count = 0;
+    catalog->tiles = NULL;
+    catalog->all_sky_root = NULL;
+    catalog->kind = STAR_CATALOG_MEMORY;
 }

src/catalog.h

@@ -9,11 +9,30 @@ typedef struct {
     double amplitude;
 } Star;
 
+enum { CATALOG_ALL_SKY_RA_TILES = 360, CATALOG_ALL_SKY_DEC_TILES = 180 };
+
 typedef struct {
     Star *stars;
     size_t count;
+    int state; /* 0: not requested, 1: loaded, -1: absent or unreadable. */
+} CatalogTile;
+
+typedef enum {
+    STAR_CATALOG_MEMORY,
+    STAR_CATALOG_ALL_SKY
+} StarCatalogKind;
+
+typedef struct {
+    Star *stars;
+    size_t count;
+    StarCatalogKind kind;
+    char *all_sky_root;
+    CatalogTile *tiles;
 } StarCatalog;
 
+typedef int (*CatalogTileVisitor)(const Star *stars, size_t count,
+                                  int fully_contained, void *context);
+
 /*
  * Synthetic lensing fixture: stars lie on the union of 10-degree longitude
  * and latitude lines, sampled every 2 degrees.  The eight longitude/hemisphere
@@ -22,6 +41,15 @@ typedef struct {
  */
 int catalog_write_octant_grid(const char *path);
 int catalog_load_csv(StarCatalog *catalog, const char *path);
+/* The directory contains tile_raRRR_decDDD.csv plus optional .done markers.
+ * Tiles are loaded only after a source triangle intersects them. */
+int catalog_load_all_sky(StarCatalog *catalog, const char *directory);
+/* Visit only 1-degree all-sky tiles that can intersect this source triangle.
+ * Call with load_missing=1 before parallel rendering, then 0 inside workers. */
+int catalog_visit_source_triangle(StarCatalog *catalog,
+                                  const double direction[3][3],
+                                  int load_missing, CatalogTileVisitor visitor,
+                                  void *context);
 void catalog_destroy(StarCatalog *catalog);
 
 #endif

src/frame.c

@@ -118,13 +118,22 @@ static double spherical_area(const double a[3], const double b[3],
                      1.0 + dot(a, b) + dot(b, c) + dot(c, a));
 }
 
+static int spherical_barycentric_weights(const double point[3], const double a[3],
+                                         const double b[3], const double c[3],
+                                         double weights[3]) {
+  const double area = spherical_area(a, b, c);
+  if (area < 1e-14)
+    return -1;
+  weights[0] = spherical_area(point, b, c) / area;
+  weights[1] = spherical_area(point, c, a) / area;
+  weights[2] = spherical_area(point, a, b) / area;
+  return 0;
+}
+
 static int spherical_barycentric(const double point[3], const double a[3],
                                  const double b[3], const double c[3],
                                  double weights[3]) {
-  const double area = spherical_area(a, b, c);
   double edge_cross[3];
-  if (area < 1e-14)
-    return -1;
   const double *corners[3] = {a, b, c};
   for (int edge = 0; edge < 3; ++edge) {
     const double *left = corners[edge];
@@ -138,10 +147,7 @@ static int spherical_barycentric(const double point[3], const double a[3],
         -1e-14 * dot(edge_cross, edge_cross))
       return -1;
   }
-  weights[0] = spherical_area(point, b, c) / area;
-  weights[1] = spherical_area(point, c, a) / area;
-  weights[2] = spherical_area(point, a, b) / area;
-  return 0;
+  return spherical_barycentric_weights(point, a, b, c, weights);
 }
 
 static int usable_triangle(const FrameLensMesh *mesh,
@@ -171,8 +177,57 @@ static int owns_source_boundary(const LensTriangle *triangle,
   return 1;
 }
 
+typedef struct {
+  const LensVertex *vertex[3];
+  const LensTriangle *triangle;
+  double *hdr;
+  int width, height;
+  double exposure, magnification;
+  const PointSpreadFunction *psf;
+  size_t images;
+} TriangleSplatContext;
+
+static int splat_catalog_tile(const Star *stars, size_t count,
+                              int fully_contained, void *opaque) {
+  TriangleSplatContext *context = opaque;
+  for (size_t s = 0; s < count; ++s) {
+    const Star *star = &stars[s];
+    double weights[3];
+    if (!fully_contained &&
+        spherical_barycentric(star->direction, context->vertex[0]->n_infinity,
+                              context->vertex[1]->n_infinity,
+                              context->vertex[2]->n_infinity, weights))
+      continue;
+    if (fully_contained) {
+      /* Only inverse-map weights remain: no per-star containment test. */
+      if (spherical_barycentric_weights(star->direction,
+                                        context->vertex[0]->n_infinity,
+                                        context->vertex[1]->n_infinity,
+                                        context->vertex[2]->n_infinity, weights))
+        return -1;
+    }
+    if (!owns_source_boundary(context->triangle, weights))
+      continue;
+    const double image_x = weights[0] * context->vertex[0]->image_x +
+                           weights[1] * context->vertex[1]->image_x +
+                           weights[2] * context->vertex[2]->image_x;
+    const double image_y = weights[0] * context->vertex[0]->image_y +
+                           weights[1] * context->vertex[1]->image_y +
+                           weights[2] * context->vertex[2]->image_y;
+    const double log_g = weights[0] * context->vertex[0]->log_frequency_ratio +
+                         weights[1] * context->vertex[1]->log_frequency_ratio +
+                         weights[2] * context->vertex[2]->log_frequency_ratio;
+    const LinearRgb color = blackbody_to_linear_rgb(star->temperature_K * exp(log_g));
+    splat_moffat(context->hdr, context->width, context->height, image_x, image_y,
+                 color, context->exposure * star->amplitude * context->magnification,
+                 context->psf);
+    ++context->images;
+  }
+  return 0;
+}
+
 static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
-                                      const StarCatalog *catalog, double *hdr,
+                                      StarCatalog *catalog, double *hdr,
                                       int width, int height, double exposure,
                                       const PointSpreadFunction *psf,
                                       size_t first_triangle,
@@ -188,42 +243,60 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
         spherical_area(vertex[0]->camera_direction, vertex[1]->camera_direction,
                        vertex[2]->camera_direction);
     const double magnification = image_area / source_area;
-    for (size_t s = 0; s < catalog->count; ++s) {
-      const Star *star = &catalog->stars[s];
-      double weights[3];
-      if (spherical_barycentric(star->direction, vertex[0]->n_infinity,
-                                vertex[1]->n_infinity, vertex[2]->n_infinity,
-                                weights))
-        continue;
-      if (!owns_source_boundary(&mesh->triangles[t], weights))
-        continue;
-      const double image_x = weights[0] * vertex[0]->image_x +
-                             weights[1] * vertex[1]->image_x +
-                             weights[2] * vertex[2]->image_x;
-      const double image_y = weights[0] * vertex[0]->image_y +
-                             weights[1] * vertex[1]->image_y +
-                             weights[2] * vertex[2]->image_y;
-      const double log_g = weights[0] * vertex[0]->log_frequency_ratio +
-                           weights[1] * vertex[1]->log_frequency_ratio +
-                           weights[2] * vertex[2]->log_frequency_ratio;
-      const LinearRgb color =
-          blackbody_to_linear_rgb(star->temperature_K * exp(log_g));
-      splat_moffat(hdr, width, height, image_x, image_y, color,
-                   exposure * star->amplitude * magnification, psf);
-      ++images;
-    }
+    const double direction[3][3] = {
+        {vertex[0]->n_infinity[0], vertex[0]->n_infinity[1], vertex[0]->n_infinity[2]},
+        {vertex[1]->n_infinity[0], vertex[1]->n_infinity[1], vertex[1]->n_infinity[2]},
+        {vertex[2]->n_infinity[0], vertex[2]->n_infinity[1], vertex[2]->n_infinity[2]}};
+    TriangleSplatContext context = {.vertex = {vertex[0], vertex[1], vertex[2]},
+                                    .triangle = &mesh->triangles[t], .hdr = hdr,
+                                    .width = width, .height = height,
+                                    .exposure = exposure, .magnification = magnification,
+                                    .psf = psf};
+    if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile,
+                                      &context) == 0)
+      images += context.images;
   }
   return images;
 }
 
+static int prefetch_catalog_tile(const Star *stars, size_t count,
+                                 int fully_contained, void *context) {
+  (void)stars;
+  (void)count;
+  (void)fully_contained;
+  (void)context;
+  return 0;
+}
+
+static void prefetch_catalog_for_mesh(const FrameLensMesh *mesh,
+                                      StarCatalog *catalog) {
+  if (catalog->kind != STAR_CATALOG_ALL_SKY)
+    return;
+  for (size_t t = 0; t < mesh->triangle_count; ++t) {
+    const LensVertex *vertex[3];
+    if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
+      continue;
+    const double direction[3][3] = {
+        {vertex[0]->n_infinity[0], vertex[0]->n_infinity[1], vertex[0]->n_infinity[2]},
+        {vertex[1]->n_infinity[0], vertex[1]->n_infinity[1], vertex[1]->n_infinity[2]},
+        {vertex[2]->n_infinity[0], vertex[2]->n_infinity[1], vertex[2]->n_infinity[2]}};
+    (void)catalog_visit_source_triangle(catalog, direction, 1,
+                                        prefetch_catalog_tile, NULL);
+  }
+}
+
 size_t frame_splat_catalog(const FrameLensMesh *mesh,
-                           const StarCatalog *catalog, double *hdr, int width,
+                           StarCatalog *catalog, double *hdr, int width,
                            int height, double exposure,
                            const PointSpreadFunction *psf) {
   if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 ||
       psf == NULL || width <= 0 || height <= 0)
     return 0;
 
+  /* Tile I/O is deliberately serial and complete before OpenMP workers start.
+   * The parallel splat pass then reads an immutable tile cache. */
+  prefetch_catalog_for_mesh(mesh, catalog);
+
   const size_t pixel_count = (size_t)width * height * 3;
   if (pixel_count > SIZE_MAX / sizeof(double) ||
       pixel_count * sizeof(double) > FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / 2)

src/frame.h

@@ -36,7 +36,7 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime,
 /* Each locally invertible escaped triangle contributes one image per contained
  * star. */
 size_t frame_splat_catalog(const FrameLensMesh *mesh,
-                           const StarCatalog *catalog, double *hdr, int width,
+                           StarCatalog *catalog, double *hdr, int width,
                            int height, double exposure,
                            const PointSpreadFunction *psf);
 void frame_draw_mesh(const FrameLensMesh *mesh, double *hdr, int width,

src/main.c

@@ -22,6 +22,7 @@ typedef struct {
   double observer_inward_speed;
   PointSpreadFunction psf;
   const char *catalog_path;
+  const char *all_sky_catalog_path;
   const char *output_path;
   const char *observer_track_path;
   const char *frames_dir;
@@ -112,6 +113,8 @@ static int parse_args(int argc, char **argv, Settings *s,
   for (int i = 1; i < argc; ++i) {
     if (!strcmp(argv[i], "--catalog") && i + 1 < argc)
       s->catalog_path = argv[++i];
+    else if (!strcmp(argv[i], "--all-sky-catalog") && i + 1 < argc)
+      s->all_sky_catalog_path = argv[++i];
     else if (!strcmp(argv[i], "--output") && i + 1 < argc)
       s->output_path = argv[++i];
     else if (!strcmp(argv[i], "--width") && i + 1 < argc &&
@@ -185,7 +188,7 @@ static int default_observer(const Settings *s, ObserverState *observer) {
 #endif
 }
 
-static int render_observer_frame(const Settings *s, const StarCatalog *catalog,
+static int render_observer_frame(const Settings *s, StarCatalog *catalog,
                                  const SpacetimeSource *spacetime,
                                  const ObserverState *observer,
                                  const char *output_path) {
@@ -214,7 +217,7 @@ static int render_observer_frame(const Settings *s, const StarCatalog *catalog,
   return result;
 }
 
-static int render_frame(const Settings *s, const StarCatalog *catalog,
+static int render_frame(const Settings *s, StarCatalog *catalog,
                         const SpacetimeSource *spacetime) {
   ObserverState observer;
   return default_observer(s, &observer) ? -1
@@ -234,7 +237,7 @@ static int frame_output_path(char path[PATH_MAX], const Settings *s,
   return written < 0 || written >= PATH_MAX ? -1 : 0;
 }
 
-static int render_movie(const Settings *s, const StarCatalog *catalog,
+static int render_movie(const Settings *s, StarCatalog *catalog,
                         const SpacetimeSource *spacetime) {
   ObserverTrack track = {0};
   Movie movie = {0};
@@ -323,7 +326,7 @@ int main(int argc, char **argv) {
   const char *write_path;
   if (parse_args(argc, argv, &settings, &write_path)) {
     fprintf(stderr,
-            "Usage: %s [--catalog PATH] [--output PATH] [--width N] [--height "
+            "Usage: %s [--catalog PATH | --all-sky-catalog DIR] [--output PATH] [--width N] [--height "
             "N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] "
             "[--exposure E] [--observer-inward-speed V] "
             "[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
@@ -341,8 +344,13 @@ int main(int argc, char **argv) {
     return write_minkowski_accel_track(&settings) == 0
                ? 0
                : (perror(settings.write_minkowski_accel_track_path), 1);
-  StarCatalog catalog;
-  if (catalog_load_csv(&catalog, settings.catalog_path)) {
+  StarCatalog catalog = {0};
+  if (settings.all_sky_catalog_path != NULL) {
+    if (catalog_load_all_sky(&catalog, settings.all_sky_catalog_path)) {
+      perror(settings.all_sky_catalog_path);
+      return 1;
+    }
+  } else if (catalog_load_csv(&catalog, settings.catalog_path)) {
     if (catalog_write_octant_grid(settings.catalog_path) ||
         catalog_load_csv(&catalog, settings.catalog_path)) {
       perror(settings.catalog_path);