diff --git a/.gitignore b/.gitignore index 6240d7d..283b0ef 100644 --- a/.gitignore +++ b/.gitignore @@ -16,6 +16,7 @@ dist/ .DS_Store .vscode/ .idea/ +.copilot-project-context.md # Generated outputs *.ppm diff --git a/Makefile b/Makefile index d54c76f..624ba14 100644 --- a/Makefile +++ b/Makefile @@ -3,13 +3,12 @@ # Unstructured data visualization tool CC = gcc -CFLAGS = -Wall -Wextra -O2 -g +CFLAGS = -Wall -Wextra -O2 -g -fopenmp # --enable-new-dtags is Linux-only, skip on macOS (Darwin) UNAME_S := $(shell uname -s) -ifneq ($(UNAME_S),Darwin) -LDFLAGS = -Wl,--enable-new-dtags -else -LDFLAGS = +LDFLAGS = -fopenmp +ifeq ($(UNAME_S),Darwin) +LDFLAGS += -Wl,--enable-new-dtags endif # NetCDF - try DKRZ system installation first, fall back to system nc-config @@ -71,7 +70,7 @@ endif BASE_CFLAGS = $(CFLAGS) $(NETCDF_CFLAGS) X11_FULL_CFLAGS = $(X11_CFLAGS) $(BASE_CFLAGS) -COMMON_LIBS = $(NETCDF_LIBS) $(NETCDF_RPATH) -lm +COMMON_LIBS = $(NETCDF_LIBS) $(NETCDF_RPATH) -lm -lpthread USHOW_LIBS = $(COMMON_LIBS) $(X11_LIBS) $(X11_RPATH) UTERM_LIBS = $(COMMON_LIBS) @@ -179,7 +178,8 @@ SRCDIR = src OBJDIR = obj BINDIR = . -COMMON_SRCS = $(SRCDIR)/kdtree.c \ +COMMON_SRCS = $(SRCDIR)/common.c \ + $(SRCDIR)/kdtree.c \ $(SRCDIR)/mesh.c \ $(SRCDIR)/regrid.c \ $(SRCDIR)/projection.c \ @@ -263,29 +263,29 @@ clean: rm -f uterm # Dependencies -$(OBJDIR)/ushow.o: $(SRCDIR)/ushow.c $(SRCDIR)/ushow.defines.h $(SRCDIR)/mesh.h \ +$(OBJDIR)/ushow.o: $(SRCDIR)/ushow.c $(SRCDIR)/us_types.h $(SRCDIR)/mesh.h \ $(SRCDIR)/regrid.h $(SRCDIR)/file_netcdf.h $(SRCDIR)/colormaps.h \ $(SRCDIR)/view.h $(SRCDIR)/interface/x_interface.h -$(OBJDIR)/uterm.o: $(SRCDIR)/uterm.c $(SRCDIR)/ushow.defines.h $(SRCDIR)/mesh.h \ +$(OBJDIR)/uterm.o: $(SRCDIR)/uterm.c $(SRCDIR)/us_types.h $(SRCDIR)/mesh.h \ $(SRCDIR)/regrid.h $(SRCDIR)/file_netcdf.h $(SRCDIR)/colormaps.h \ $(SRCDIR)/term_render_mode.h \ $(SRCDIR)/view.h $(OBJDIR)/term_render_mode.o: $(SRCDIR)/term_render_mode.c $(SRCDIR)/term_render_mode.h $(OBJDIR)/kdtree.o: $(SRCDIR)/kdtree.c $(SRCDIR)/kdtree.h -$(OBJDIR)/mesh.o: $(SRCDIR)/mesh.c $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h +$(OBJDIR)/mesh.o: $(SRCDIR)/mesh.c $(SRCDIR)/mesh.h $(SRCDIR)/us_types.h $(OBJDIR)/regrid.o: $(SRCDIR)/regrid.c $(SRCDIR)/regrid.h $(SRCDIR)/mesh.h \ - $(SRCDIR)/kdtree.h $(SRCDIR)/projection.h $(SRCDIR)/ushow.defines.h -$(OBJDIR)/projection.o: $(SRCDIR)/projection.c $(SRCDIR)/projection.h $(SRCDIR)/ushow.defines.h -$(OBJDIR)/file_netcdf.o: $(SRCDIR)/file_netcdf.c $(SRCDIR)/file_netcdf.h $(SRCDIR)/ushow.defines.h -$(OBJDIR)/file_mitgcm.o: $(SRCDIR)/file_mitgcm.c $(SRCDIR)/file_mitgcm.h $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h -$(OBJDIR)/colormaps.o: $(SRCDIR)/colormaps.c $(SRCDIR)/colormaps.h $(SRCDIR)/ushow.defines.h + $(SRCDIR)/kdtree.h $(SRCDIR)/projection.h $(SRCDIR)/us_types.h +$(OBJDIR)/projection.o: $(SRCDIR)/projection.c $(SRCDIR)/projection.h $(SRCDIR)/us_types.h +$(OBJDIR)/file_netcdf.o: $(SRCDIR)/file_netcdf.c $(SRCDIR)/file_netcdf.h $(SRCDIR)/us_types.h +$(OBJDIR)/file_mitgcm.o: $(SRCDIR)/file_mitgcm.c $(SRCDIR)/file_mitgcm.h $(SRCDIR)/mesh.h $(SRCDIR)/us_types.h +$(OBJDIR)/colormaps.o: $(SRCDIR)/colormaps.c $(SRCDIR)/colormaps.h $(SRCDIR)/us_types.h $(OBJDIR)/view.o: $(SRCDIR)/view.c $(SRCDIR)/view.h $(SRCDIR)/file_netcdf.h \ - $(SRCDIR)/regrid.h $(SRCDIR)/colormaps.h $(SRCDIR)/ushow.defines.h + $(SRCDIR)/regrid.h $(SRCDIR)/colormaps.h $(SRCDIR)/us_types.h $(OBJDIR)/interface/x_interface.o: $(SRCDIR)/interface/x_interface.c \ $(SRCDIR)/interface/x_interface.h \ $(SRCDIR)/interface/colorbar.h \ $(SRCDIR)/interface/timeseries_popup.h \ - $(SRCDIR)/ushow.defines.h + $(SRCDIR)/us_types.h $(OBJDIR)/interface/colorbar.o: $(SRCDIR)/interface/colorbar.c \ $(SRCDIR)/interface/colorbar.h $(SRCDIR)/colormaps.h $(OBJDIR)/interface/range_popup.o: $(SRCDIR)/interface/range_popup.c \ @@ -295,25 +295,25 @@ $(OBJDIR)/interface/range_utils.o: $(SRCDIR)/interface/range_utils.c \ $(SRCDIR)/interface/range_utils.h $(OBJDIR)/interface/timeseries_popup.o: $(SRCDIR)/interface/timeseries_popup.c \ $(SRCDIR)/interface/timeseries_popup.h \ - $(SRCDIR)/ushow.defines.h + $(SRCDIR)/us_types.h # Zarr dependencies (when WITH_ZARR is set) ifdef WITH_ZARR -$(OBJDIR)/file_zarr.o: $(SRCDIR)/file_zarr.c $(SRCDIR)/file_zarr.h $(SRCDIR)/ushow.defines.h \ +$(OBJDIR)/file_zarr.o: $(SRCDIR)/file_zarr.c $(SRCDIR)/file_zarr.h $(SRCDIR)/us_types.h \ $(SRCDIR)/cJSON/cJSON.h $(OBJDIR)/cJSON/cJSON.o: $(SRCDIR)/cJSON/cJSON.c $(SRCDIR)/cJSON/cJSON.h endif # GRIB dependencies (when WITH_GRIB is set) ifdef WITH_GRIB -$(OBJDIR)/file_grib.o: $(SRCDIR)/file_grib.c $(SRCDIR)/file_grib.h $(SRCDIR)/ushow.defines.h +$(OBJDIR)/file_grib.o: $(SRCDIR)/file_grib.c $(SRCDIR)/file_grib.h $(SRCDIR)/us_types.h endif # YAC dependencies (when WITH_YAC is set) ifdef WITH_YAC $(OBJDIR)/regrid_yac.o: $(SRCDIR)/regrid_yac.c $(SRCDIR)/regrid_yac.h \ - $(SRCDIR)/healpix.h $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h -$(OBJDIR)/healpix.o: $(SRCDIR)/healpix.c $(SRCDIR)/healpix.h $(SRCDIR)/ushow.defines.h + $(SRCDIR)/healpix.h $(SRCDIR)/mesh.h $(SRCDIR)/us_types.h +$(OBJDIR)/healpix.o: $(SRCDIR)/healpix.c $(SRCDIR)/healpix.h $(SRCDIR)/us_types.h endif # Print configuration diff --git a/src/colormaps.h b/src/colormaps.h index 07ca5c9..9bceacd 100644 --- a/src/colormaps.h +++ b/src/colormaps.h @@ -5,7 +5,7 @@ #ifndef COLORMAPS_H #define COLORMAPS_H -#include "ushow.defines.h" +#include "us_types.h" /* * Initialize colormaps (load built-in colormaps). diff --git a/src/common.c b/src/common.c new file mode 100644 index 0000000..75bc4da --- /dev/null +++ b/src/common.c @@ -0,0 +1,134 @@ +/* + * common.c - Shared utilities for ushow and uterm + */ + +#include "common.h" +#include "kdtree.h" + +#include +#include +#include +#include +#ifdef _OPENMP +#include +#endif + +#define DEFAULT_THREADS 4 + +void apply_thread_settings(int user_threads) { + if (user_threads > 0) { + kdtree_set_max_threads(user_threads); +#ifdef _OPENMP + omp_set_num_threads(user_threads); +#endif + } else if (!getenv("OMP_NUM_THREADS")) { +#ifdef _OPENMP + omp_set_num_threads(DEFAULT_THREADS); +#endif + } +} + +int parse_time_units(const char *units, double *unit_seconds, + int *y, int *mo, int *d, int *h, int *mi, double *sec) { + if (!units || !unit_seconds || !y || !mo || !d || !h || !mi || !sec) return 0; + + const char *since = strstr(units, "since"); + if (!since) return 0; + + char unit_buf[32] = {0}; + if (sscanf(units, "%31s", unit_buf) != 1) return 0; + + /* Normalize unit token to lower case */ + for (char *p = unit_buf; *p; ++p) { + if (*p >= 'A' && *p <= 'Z') *p = (char)(*p - 'A' + 'a'); + } + + if (strcmp(unit_buf, "seconds") == 0 || strcmp(unit_buf, "second") == 0 || + strcmp(unit_buf, "secs") == 0 || strcmp(unit_buf, "sec") == 0 || strcmp(unit_buf, "s") == 0) { + *unit_seconds = 1.0; + } else if (strcmp(unit_buf, "minutes") == 0 || strcmp(unit_buf, "minute") == 0 || + strcmp(unit_buf, "mins") == 0 || strcmp(unit_buf, "min") == 0) { + *unit_seconds = 60.0; + } else if (strcmp(unit_buf, "hours") == 0 || strcmp(unit_buf, "hour") == 0 || + strcmp(unit_buf, "hrs") == 0 || strcmp(unit_buf, "hr") == 0) { + *unit_seconds = 3600.0; + } else if (strcmp(unit_buf, "days") == 0 || strcmp(unit_buf, "day") == 0) { + *unit_seconds = 86400.0; + } else { + return 0; + } + + /* Parse origin date/time after "since" */ + const char *p = since + 5; + while (*p == ' ') p++; + int n = sscanf(p, "%d-%d-%d %d:%d:%lf", y, mo, d, h, mi, sec); + if (n < 3) return 0; + if (n == 3) { *h = 0; *mi = 0; *sec = 0.0; } + return 1; +} + +int64_t days_from_civil(int y, unsigned m, unsigned d) { + y -= m <= 2; + const int era = (y >= 0 ? y : y - 399) / 400; + const unsigned yoe = (unsigned)(y - era * 400); + const unsigned doy = (153 * (m + (m > 2 ? -3 : 9)) + 2) / 5 + d - 1; + const unsigned doe = yoe * 365 + yoe / 4 - yoe / 100 + doy; + return (int64_t)(era * 146097 + (int)doe - 719468); +} + +void civil_from_days(int64_t z, int *y, unsigned *m, unsigned *d) { + z += 719468; + const int era = (z >= 0 ? z : z - 146096) / 146097; + const unsigned doe = (unsigned)(z - era * 146097); + const unsigned yoe = (doe - doe / 1460 + doe / 36524 - doe / 146096) / 365; + int y_tmp = (int)(yoe) + era * 400; + const unsigned doy = doe - (365 * yoe + yoe / 4 - yoe / 100); + const unsigned mp = (5 * doy + 2) / 153; + const unsigned d_tmp = doy - (153 * mp + 2) / 5 + 1; + const unsigned m_tmp = mp + (mp < 10 ? 3 : -9); + y_tmp += (m_tmp <= 2); + + *y = y_tmp; + *m = m_tmp; + *d = d_tmp; +} + +int format_time_from_units(char *out, size_t outlen, double value, const char *units) { + if (!out || outlen == 0) return 0; + + double unit_seconds = 0.0; + int y, mo, d, h, mi; + double sec; + if (!parse_time_units(units, &unit_seconds, &y, &mo, &d, &h, &mi, &sec)) { + return 0; + } + + int64_t days = days_from_civil(y, (unsigned)mo, (unsigned)d); + double total_sec = (double)days * 86400.0 + (double)h * 3600.0 + (double)mi * 60.0 + sec; + total_sec += value * unit_seconds; + + int64_t out_days = (int64_t)(total_sec / 86400.0); + double rem = total_sec - (double)out_days * 86400.0; + if (rem < 0) { + rem += 86400.0; + out_days -= 1; + } + + int out_y; + unsigned out_m, out_d; + civil_from_days(out_days, &out_y, &out_m, &out_d); + + int out_h = (int)(rem / 3600.0); + rem -= out_h * 3600.0; + int out_mi = (int)(rem / 60.0); + rem -= out_mi * 60.0; + int out_s = (int)(rem + 0.5); + + if (out_h == 0 && out_mi == 0 && out_s == 0) { + snprintf(out, outlen, "%04d-%02u-%02u", out_y, out_m, out_d); + } else { + snprintf(out, outlen, "%04d-%02u-%02u %02d:%02d:%02d", + out_y, out_m, out_d, out_h, out_mi, out_s); + } + return 1; +} diff --git a/src/common.h b/src/common.h new file mode 100644 index 0000000..8c0e3cc --- /dev/null +++ b/src/common.h @@ -0,0 +1,41 @@ +/* + * common.h - Shared utilities for ushow and uterm + */ + +#ifndef COMMON_H +#define COMMON_H + +#include +#include + +/* + * Apply thread settings for kdtree and OpenMP. + * user_threads > 0: use that value for both kdtree and omp_set_num_threads. + * user_threads <= 0: leave kdtree default; if OMP_NUM_THREADS is unset, + * set omp_set_num_threads to 4. + */ +void apply_thread_settings(int user_threads); + +/* + * Parse a CF-convention "units since origin" string. + * Returns 1 on success, 0 on failure. + * unit_seconds: conversion factor (e.g. 86400 for days) + * y,mo,d,h,mi,sec: parsed origin date/time components + */ +int parse_time_units(const char *units, double *unit_seconds, + int *y, int *mo, int *d, int *h, int *mi, double *sec); + +/* + * Calendar arithmetic helpers (proleptic Gregorian). + */ +int64_t days_from_civil(int y, unsigned m, unsigned d); +void civil_from_days(int64_t z, int *y, unsigned *m, unsigned *d); + +/* + * Format a CF time value as "YYYY-MM-DD HH:MM:SS" (or "YYYY-MM-DD" when + * the time-of-day is exactly midnight). + * Returns 1 on success, 0 on failure. + */ +int format_time_from_units(char *out, size_t outlen, double value, const char *units); + +#endif /* COMMON_H */ diff --git a/src/file_grib.h b/src/file_grib.h index c733b90..b6cb08f 100644 --- a/src/file_grib.h +++ b/src/file_grib.h @@ -8,7 +8,7 @@ #ifndef FILE_GRIB_H #define FILE_GRIB_H -#include "ushow.defines.h" +#include "us_types.h" #ifdef HAVE_GRIB diff --git a/src/file_mitgcm.c b/src/file_mitgcm.c index e4f8827..4076732 100644 --- a/src/file_mitgcm.c +++ b/src/file_mitgcm.c @@ -407,7 +407,7 @@ USFile *mitgcm_open(const char *path) { /* Allocate store */ MitgcmStore *store = calloc(1, sizeof(MitgcmStore)); if (!store) return NULL; - strncpy(store->directory, directory, MAX_NAME_LEN - 1); + snprintf(store->directory, sizeof(store->directory), "%s", directory); store->nx = nx; store->ny = ny; store->nz = nz; @@ -473,7 +473,7 @@ USFile *mitgcm_open(const char *path) { free(store); return NULL; } f->file_type = FILE_TYPE_MITGCM; - strncpy(f->filename, directory, MAX_NAME_LEN - 1); + snprintf(f->filename, sizeof(f->filename), "%s", directory); f->mitgcm_data = store; return f; @@ -487,7 +487,7 @@ USMesh *mitgcm_create_mesh(USFile *file) { size_t n_points = (size_t)ny * nx; /* Read XC.data and YC.data */ - char xc_path[MAX_NAME_LEN], yc_path[MAX_NAME_LEN]; + char xc_path[MAX_NAME_LEN + 16], yc_path[MAX_NAME_LEN + 16]; snprintf(xc_path, sizeof(xc_path), "%s/XC.data", store->directory); snprintf(yc_path, sizeof(yc_path), "%s/YC.data", store->directory); @@ -637,7 +637,8 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) { if (strcmp(prefixes[i], prefix) == 0) { found = 1; break; } } if (!found && n_prefixes < 64) { - strncpy(prefixes[n_prefixes++], prefix, MAX_NAME_LEN - 1); + snprintf(prefixes[n_prefixes], MAX_NAME_LEN, "%s", prefix); + n_prefixes++; } } closedir(dir); @@ -659,7 +660,7 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) { if (!iterations || n_iters == 0) { free(iterations); continue; } /* Parse meta from first iteration to get field info */ - char meta_path[MAX_NAME_LEN]; + char meta_path[MAX_NAME_LEN * 2 + 32]; snprintf(meta_path, sizeof(meta_path), "%s/%s.%010d.meta", store->directory, prefixes[pi], iterations[0]); @@ -683,9 +684,9 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) { /* Name */ if (minfo.n_fields > 0 && minfo.fields[fi][0]) { - strncpy(var->name, minfo.fields[fi], MAX_NAME_LEN - 1); + snprintf(var->name, sizeof(var->name), "%s", minfo.fields[fi]); } else { - strncpy(var->name, prefixes[pi], MAX_NAME_LEN - 1); + snprintf(var->name, sizeof(var->name), "%s", prefixes[pi]); } var->mesh = m; @@ -726,7 +727,7 @@ USVar *mitgcm_scan_variables(USFile *f, USMesh *m) { /* Private data */ MitgcmVarData *vd = calloc(1, sizeof(MitgcmVarData)); if (!vd) { free(var); continue; } - strncpy(vd->prefix, prefixes[pi], MAX_NAME_LEN - 1); + snprintf(vd->prefix, sizeof(vd->prefix), "%s", prefixes[pi]); vd->field_index = fi; vd->iterations = malloc(n_iters * sizeof(int)); if (vd->iterations) { @@ -801,7 +802,7 @@ int mitgcm_read_slice(USVar *var, size_t time_idx, size_t depth_idx, float *data if ((int)time_idx >= vd->n_iterations) return -1; /* Build data file path */ - char data_path[MAX_NAME_LEN]; + char data_path[MAX_NAME_LEN * 2 + 32]; snprintf(data_path, sizeof(data_path), "%s/%s.%010d.data", store->directory, vd->prefix, vd->iterations[time_idx]); @@ -1054,7 +1055,7 @@ int mitgcm_read_timeseries(USVar *var, size_t node_idx, size_t depth_idx, times[t] = (double)vd->iterations[t]; /* Build path and read single value */ - char data_path[MAX_NAME_LEN]; + char data_path[MAX_NAME_LEN * 2 + 32]; snprintf(data_path, sizeof(data_path), "%s/%s.%010d.data", store->directory, vd->prefix, vd->iterations[t]); diff --git a/src/file_mitgcm.h b/src/file_mitgcm.h index 14d1bf7..01e8d97 100644 --- a/src/file_mitgcm.h +++ b/src/file_mitgcm.h @@ -8,7 +8,7 @@ #ifndef FILE_MITGCM_H #define FILE_MITGCM_H -#include "ushow.defines.h" +#include "us_types.h" /* * Check if a path looks like MITgcm data (directory with .meta files, diff --git a/src/file_netcdf.h b/src/file_netcdf.h index c9fd082..dbbd03a 100644 --- a/src/file_netcdf.h +++ b/src/file_netcdf.h @@ -5,7 +5,7 @@ #ifndef FILE_NETCDF_H #define FILE_NETCDF_H -#include "ushow.defines.h" +#include "us_types.h" /* * Open a NetCDF file and create file structure. diff --git a/src/file_zarr.h b/src/file_zarr.h index e5ad073..4f4d037 100644 --- a/src/file_zarr.h +++ b/src/file_zarr.h @@ -8,7 +8,7 @@ #ifndef FILE_ZARR_H #define FILE_ZARR_H -#include "ushow.defines.h" +#include "us_types.h" #ifdef HAVE_ZARR diff --git a/src/healpix.h b/src/healpix.h index 74f382d..c3a80b8 100644 --- a/src/healpix.h +++ b/src/healpix.h @@ -9,7 +9,7 @@ #ifndef HEALPIX_H #define HEALPIX_H -#include "ushow.defines.h" +#include "us_types.h" /* * Generate triangular element connectivity for a mesh whose points are diff --git a/src/interface/timeseries_popup.h b/src/interface/timeseries_popup.h index 9a65b66..2ea98d4 100644 --- a/src/interface/timeseries_popup.h +++ b/src/interface/timeseries_popup.h @@ -9,7 +9,7 @@ #define TIMESERIES_POPUP_H #include -#include "../ushow.defines.h" +#include "../us_types.h" /* * Initialize the timeseries popup widgets. diff --git a/src/interface/x_interface.h b/src/interface/x_interface.h index a0dd6ba..c96f107 100644 --- a/src/interface/x_interface.h +++ b/src/interface/x_interface.h @@ -5,7 +5,7 @@ #ifndef X_INTERFACE_H #define X_INTERFACE_H -#include "../ushow.defines.h" +#include "../us_types.h" /* * Set light theme (call before x_init). Default is dark theme. diff --git a/src/kdtree.c b/src/kdtree.c index 488f6c0..27efda8 100644 --- a/src/kdtree.c +++ b/src/kdtree.c @@ -1,22 +1,61 @@ /* - * kdtree.c - Simple KDTree implementation for 3D nearest-neighbor queries + * kdtree.c - Advanced KDTree implementation for 3D nearest-neighbor queries * - * Uses median-split construction and recursive nearest-neighbor search. - * Optimized for the case of building once and querying many times. + * Inspired by the algorithms used in libkdtree by Jörg Dietrich + * (https://github.com/joergdietrich/libkdtree, LGPL-2.0). + * No code was copied; this is an independent reimplementation using + * standard kd-tree techniques (hyperrectangle-bounded nodes, + * parallel construction, multi-dimensional pruning) adapted for + * 3D double-precision coordinates with the existing public API. */ +#define _GNU_SOURCE /* for qsort_r on Linux */ #include "kdtree.h" #include #include #include #include +#include +#include #define KDTREE_DIM 3 -/* KDTree node */ +/* Minimum number of points before spawning a new thread */ +#define KDTREE_THREAD_THRESHOLD 4096 + +/* Default thread count when neither set by user nor OMP_NUM_THREADS */ +#define KDTREE_DEFAULT_THREADS 4 + +/* Module-level thread limit: -1 means "not explicitly set" */ +static int g_max_threads = -1; + +void kdtree_set_max_threads(int n) { + if (n > 0) g_max_threads = n; +} + +/* Resolve effective thread count: + * 1. Explicit call to kdtree_set_max_threads() + * 2. OMP_NUM_THREADS environment variable + * 3. Number of online CPUs (sysconf) + * 4. KDTREE_DEFAULT_THREADS + */ +static int get_max_threads(void) { + if (g_max_threads > 0) return g_max_threads; + const char *env = getenv("OMP_NUM_THREADS"); + if (env) { + int n = atoi(env); + if (n > 0) return n; + } + return KDTREE_DEFAULT_THREADS; +} + +/* KDTree node with hyperrectangle bounds (from libkdtree design) */ typedef struct KDNode { - size_t idx; /* Original point index */ - double point[KDTREE_DIM]; + size_t idx; /* Original point index */ + double location[KDTREE_DIM]; /* Node position */ + double hr_min[KDTREE_DIM]; /* Hyperrectangle min bounds */ + double hr_max[KDTREE_DIM]; /* Hyperrectangle max bounds */ + int split; /* Split axis */ struct KDNode *left; struct KDNode *right; } KDNode; @@ -28,56 +67,166 @@ struct KDTree { double *points; /* Copy of original points */ }; -/* Comparison context for qsort */ +/* Thread argument for parallel tree construction (from libkdtree design) */ +typedef struct { + const double *points; + size_t *indices; + size_t n; + double hr_min[KDTREE_DIM]; + double hr_max[KDTREE_DIM]; + int depth; + int max_threads; +} KDBuildArg; + +/* ------------------------------------------------------------------- + * Comparison for qsort_r (thread-safe, no global state) + * ------------------------------------------------------------------- */ + typedef struct { const double *points; int axis; -} SortContext; +} SortCtx; -static SortContext sort_ctx; - -/* Comparison function for sorting indices by coordinate */ -static int compare_by_axis(const void *a, const void *b) { +static int compare_by_axis_r(const void *a, const void *b, void *arg) { + SortCtx *ctx = (SortCtx *)arg; size_t ia = *(const size_t *)a; size_t ib = *(const size_t *)b; - double va = sort_ctx.points[ia * KDTREE_DIM + sort_ctx.axis]; - double vb = sort_ctx.points[ib * KDTREE_DIM + sort_ctx.axis]; + double va = ctx->points[ia * KDTREE_DIM + ctx->axis]; + double vb = ctx->points[ib * KDTREE_DIM + ctx->axis]; if (va < vb) return -1; - if (va > vb) return 1; + if (va > vb) return 1; return 0; } -/* Build KDTree recursively */ -static KDNode *build_tree(const double *points, size_t *indices, - size_t n, int depth) { +/* ------------------------------------------------------------------- + * Squared Euclidean distance + * ------------------------------------------------------------------- */ + +static inline double dist_sq(const double *a, const double *b) { + double dx = a[0] - b[0]; + double dy = a[1] - b[1]; + double dz = a[2] - b[2]; + return dx*dx + dy*dy + dz*dz; +} + +/* ------------------------------------------------------------------- + * Minimum squared distance from a point to a hyperrectangle. + * This is the key pruning improvement: instead of checking only the + * split-axis distance (simple impl) or the distance to the further + * node's split-axis bounds (original libkdtree), we compute the true + * minimum distance across all dimensions. + * ------------------------------------------------------------------- */ + +static inline double min_dist_sq_to_hr(const double *query, + const double *hr_min, + const double *hr_max) { + double dsq = 0.0; + for (int d = 0; d < KDTREE_DIM; d++) { + if (query[d] < hr_min[d]) { + double diff = hr_min[d] - query[d]; + dsq += diff * diff; + } else if (query[d] > hr_max[d]) { + double diff = query[d] - hr_max[d]; + dsq += diff * diff; + } + } + return dsq; +} + +/* ------------------------------------------------------------------- + * Tree construction (adapted from libkdtree: kd_doBuildTree) + * + * Each node stores the bounding hyperrectangle for its subtree, + * enabling tighter pruning during nearest-neighbor search. + * Construction can be parallelized across POSIX threads. + * ------------------------------------------------------------------- */ + +static KDNode *build_tree(const double *points, size_t *indices, size_t n, + const double *hr_min, const double *hr_max, + int depth, int max_threads); +static void *build_tree_thread(void *arg); + +static KDNode *build_tree(const double *points, size_t *indices, size_t n, + const double *hr_min, const double *hr_max, + int depth, int max_threads) { if (n == 0) return NULL; int axis = depth % KDTREE_DIM; - /* Sort indices by current axis */ - sort_ctx.points = points; - sort_ctx.axis = axis; - qsort(indices, n, sizeof(size_t), compare_by_axis); + /* Thread-safe sort by current axis (cf. libkdtree pmergesort) */ + SortCtx ctx = { .points = points, .axis = axis }; + qsort_r(indices, n, sizeof(size_t), compare_by_axis_r, &ctx); - /* Find median */ size_t median = n / 2; - /* Create node */ + /* Allocate node with hyperrectangle bounds (cf. libkdtree kd_allocNode) */ KDNode *node = malloc(sizeof(KDNode)); if (!node) return NULL; - node->idx = indices[median]; - node->point[0] = points[node->idx * KDTREE_DIM + 0]; - node->point[1] = points[node->idx * KDTREE_DIM + 1]; - node->point[2] = points[node->idx * KDTREE_DIM + 2]; - - /* Build subtrees */ - node->left = build_tree(points, indices, median, depth + 1); - node->right = build_tree(points, indices + median + 1, n - median - 1, depth + 1); + node->idx = indices[median]; + node->split = axis; + memcpy(node->location, &points[node->idx * KDTREE_DIM], + KDTREE_DIM * sizeof(double)); + memcpy(node->hr_min, hr_min, KDTREE_DIM * sizeof(double)); + memcpy(node->hr_max, hr_max, KDTREE_DIM * sizeof(double)); + node->left = NULL; + node->right = NULL; + + if (n == 1) return node; + + /* Child hyperrectangle bounds split at the median's coordinate + * (cf. libkdtree: tmpMaxLeft[sortaxis] = node->location[sortaxis]) */ + double left_max[KDTREE_DIM], right_min[KDTREE_DIM]; + memcpy(left_max, hr_max, KDTREE_DIM * sizeof(double)); + memcpy(right_min, hr_min, KDTREE_DIM * sizeof(double)); + left_max[axis] = node->location[axis]; + right_min[axis] = node->location[axis]; + + /* Parallel construction (cf. libkdtree: pthread_create in kd_doBuildTree) */ + if (max_threads > 1 && n > KDTREE_THREAD_THRESHOLD) { + KDBuildArg left_arg = { + .points = points, + .indices = indices, + .n = median, + .depth = depth + 1, + .max_threads = max_threads / 2 + }; + memcpy(left_arg.hr_min, hr_min, KDTREE_DIM * sizeof(double)); + memcpy(left_arg.hr_max, left_max, KDTREE_DIM * sizeof(double)); + + pthread_t tid; + if (pthread_create(&tid, NULL, build_tree_thread, &left_arg) == 0) { + /* Build right subtree in this thread */ + node->right = build_tree(points, indices + median + 1, + n - median - 1, right_min, hr_max, + depth + 1, max_threads / 2); + pthread_join(tid, (void **)&node->left); + } else { + /* pthread_create failed – fall back to sequential */ + goto sequential; + } + } else { +sequential: + node->left = build_tree(points, indices, median, + hr_min, left_max, depth + 1, max_threads); + node->right = build_tree(points, indices + median + 1, + n - median - 1, right_min, hr_max, + depth + 1, max_threads); + } return node; } +static void *build_tree_thread(void *arg) { + KDBuildArg *d = (KDBuildArg *)arg; + return build_tree(d->points, d->indices, d->n, + d->hr_min, d->hr_max, d->depth, d->max_threads); +} + +/* ------------------------------------------------------------------- + * Public: create tree + * ------------------------------------------------------------------- */ + KDTree *kdtree_create(const double *points, size_t n_points) { if (!points || n_points == 0) return NULL; @@ -101,12 +250,26 @@ KDTree *kdtree_create(const double *points, size_t n_points) { free(tree); return NULL; } + + /* Compute bounding hyperrectangle (cf. libkdtree: min/max args) */ + double bb_min[KDTREE_DIM], bb_max[KDTREE_DIM]; + for (int d = 0; d < KDTREE_DIM; d++) { + bb_min[d] = DBL_MAX; + bb_max[d] = -DBL_MAX; + } for (size_t i = 0; i < n_points; i++) { indices[i] = i; + for (int d = 0; d < KDTREE_DIM; d++) { + double v = points[i * KDTREE_DIM + d]; + if (v < bb_min[d]) bb_min[d] = v; + if (v > bb_max[d]) bb_max[d] = v; + } } - /* Build tree */ - tree->root = build_tree(tree->points, indices, n_points, 0); + /* Resolve thread count: CLI > OMP_NUM_THREADS > default */ + int max_threads = get_max_threads(); + tree->root = build_tree(tree->points, indices, n_points, + bb_min, bb_max, 0, max_threads); free(indices); @@ -119,39 +282,45 @@ KDTree *kdtree_create(const double *points, size_t n_points) { return tree; } -/* Squared Euclidean distance */ -static inline double dist_sq(const double *a, const double *b) { - double dx = a[0] - b[0]; - double dy = a[1] - b[1]; - double dz = a[2] - b[2]; - return dx*dx + dy*dy + dz*dz; -} +/* ------------------------------------------------------------------- + * Nearest-neighbor search (adapted from libkdtree: kd_nearest) + * + * Uses hyperrectangle bounds to prune the further subtree: if the + * minimum distance from the query to the further child's bounding + * box exceeds the current best, the entire subtree is skipped. + * ------------------------------------------------------------------- */ -/* Recursive nearest neighbor search */ static void search_nearest(const KDNode *node, const double *query, - int depth, size_t *best_idx, double *best_dist_sq) { + size_t *best_idx, double *best_dist_sq) { if (!node) return; /* Check current node */ - double d = dist_sq(node->point, query); + double d = dist_sq(node->location, query); if (d < *best_dist_sq) { *best_dist_sq = d; *best_idx = node->idx; } - /* Determine which subtree to search first */ - int axis = depth % KDTREE_DIM; - double diff = query[axis] - node->point[axis]; - - KDNode *first = (diff < 0) ? node->left : node->right; - KDNode *second = (diff < 0) ? node->right : node->left; + /* Determine nearer / further subtree (cf. libkdtree kd_nearest) */ + const KDNode *nearer, *further; + if (query[node->split] < node->location[node->split]) { + nearer = node->left; + further = node->right; + } else { + nearer = node->right; + further = node->left; + } /* Search closer subtree first */ - search_nearest(first, query, depth + 1, best_idx, best_dist_sq); - - /* Check if we need to search the other subtree */ - if (diff * diff < *best_dist_sq) { - search_nearest(second, query, depth + 1, best_idx, best_dist_sq); + search_nearest(nearer, query, best_idx, best_dist_sq); + + /* Prune using minimum distance to further child's hyperrectangle. + * This is tighter than the single-axis check in the simple + * implementation, and uses all dimensions unlike the original + * libkdtree which only checked the split axis. */ + if (further && min_dist_sq_to_hr(query, further->hr_min, + further->hr_max) < *best_dist_sq) { + search_nearest(further, query, best_idx, best_dist_sq); } } @@ -163,16 +332,19 @@ void kdtree_query_nearest(const KDTree *tree, const double *query, return; } - *nn_dist = DBL_MAX; + double best_dist_sq = DBL_MAX; *nn_idx = 0; - search_nearest(tree->root, query, 0, nn_idx, nn_dist); + search_nearest(tree->root, query, nn_idx, &best_dist_sq); /* Return actual distance (not squared) */ - *nn_dist = sqrt(*nn_dist); + *nn_dist = sqrt(best_dist_sq); } -/* Free tree recursively */ +/* ------------------------------------------------------------------- + * Cleanup (cf. libkdtree: kd_destroyTree / kd_freeNode) + * ------------------------------------------------------------------- */ + static void free_node(KDNode *node) { if (!node) return; free_node(node->left); diff --git a/src/kdtree.h b/src/kdtree.h index f91ea66..8ed6985 100644 --- a/src/kdtree.h +++ b/src/kdtree.h @@ -39,4 +39,12 @@ void kdtree_free(KDTree *tree); */ size_t kdtree_size(const KDTree *tree); +/* + * Set the maximum number of threads for parallel tree construction. + * If never called, kdtree_create uses OMP_NUM_THREADS (if set), + * otherwise falls back to KDTREE_DEFAULT_THREADS (4). + * A value <= 0 is ignored. + */ +void kdtree_set_max_threads(int n); + #endif /* KDTREE_H */ diff --git a/src/mesh.h b/src/mesh.h index 536c35c..7280128 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -5,7 +5,7 @@ #ifndef MESH_H #define MESH_H -#include "ushow.defines.h" +#include "us_types.h" /* * Convert longitude/latitude to Cartesian on unit sphere. diff --git a/src/projection.h b/src/projection.h index 5c5cc71..fecf71f 100644 --- a/src/projection.h +++ b/src/projection.h @@ -5,7 +5,7 @@ #ifndef PROJECTION_H #define PROJECTION_H -#include "ushow.defines.h" +#include "us_types.h" /* * Initialize target config with global equirectangular defaults. diff --git a/src/regrid.h b/src/regrid.h index f511f9d..d5093a3 100644 --- a/src/regrid.h +++ b/src/regrid.h @@ -5,7 +5,7 @@ #ifndef REGRID_H #define REGRID_H -#include "ushow.defines.h" +#include "us_types.h" #include "mesh.h" /* diff --git a/src/regrid_yac.h b/src/regrid_yac.h index 5800ed3..09ae364 100644 --- a/src/regrid_yac.h +++ b/src/regrid_yac.h @@ -11,7 +11,7 @@ #ifdef HAVE_YAC -#include "ushow.defines.h" +#include "us_types.h" #include "mesh.h" typedef enum { diff --git a/src/ushow.defines.h b/src/us_types.h similarity index 92% rename from src/ushow.defines.h rename to src/us_types.h index ed8afe5..596d6a5 100644 --- a/src/ushow.defines.h +++ b/src/us_types.h @@ -1,11 +1,11 @@ /* - * ushow.defines.h - Data structures and constants for ushow + * us_types.h - Shared data structures and constants for ushow/uterm * * Unstructured data visualization tool */ -#ifndef USHOW_DEFINES_H -#define USHOW_DEFINES_H +#ifndef US_TYPES_H +#define US_TYPES_H #include @@ -269,21 +269,6 @@ struct USView { int frame_delay_ms; }; -/* Global options */ -typedef struct { - int debug; - double influence_radius; /* Regrid influence radius in meters */ - double target_resolution; /* Target grid resolution in degrees */ - char mesh_file[MAX_NAME_LEN]; /* Separate mesh file path */ - int frame_delay_ms; /* Animation speed */ - int polygon_only; /* Skip regridding, polygon mode only */ - USTargetConfig target_config; /* Target grid configuration */ -#ifdef HAVE_YAC - int yac_method; /* YAC interpolation method (-1 = disabled) */ - int yac_3d; /* Fractional fill-value masking */ -#endif -} USOptions; - /* Dimension info for display */ typedef struct { char name[MAX_NAME_LEN]; /* Dimension name (e.g., "time", "depth") */ @@ -320,4 +305,4 @@ typedef struct { USColor *colors; } USColormap; -#endif /* USHOW_DEFINES_H */ +#endif /* US_TYPES_H */ diff --git a/src/ushow.c b/src/ushow.c index f301ec8..2508f4e 100644 --- a/src/ushow.c +++ b/src/ushow.c @@ -4,7 +4,7 @@ * Unstructured data visualization tool */ -#include "ushow.defines.h" +#include "us_types.h" #include "mesh.h" #include "regrid.h" #ifdef HAVE_YAC @@ -19,6 +19,7 @@ #include "file_grib.h" #endif #include "kdtree.h" +#include "common.h" #include "colormaps.h" #include "view.h" #include "interface/x_interface.h" @@ -50,6 +51,21 @@ static USDimInfo *current_dim_info = NULL; static int n_current_dims = 0; /* Options */ +typedef struct { + int debug; + double influence_radius; /* Regrid influence radius in meters */ + double target_resolution; /* Target grid resolution in degrees */ + char mesh_file[MAX_NAME_LEN]; /* Separate mesh file path */ + int frame_delay_ms; /* Animation speed */ + int polygon_only; /* Skip regridding, polygon mode only */ + USTargetConfig target_config; /* Target grid configuration */ + int user_threads; /* CLI --threads value (0 = not set) */ +#ifdef HAVE_YAC + int yac_method; /* YAC interpolation method (-1 = disabled) */ + int yac_3d; /* Fractional fill-value masking */ +#endif +} USOptions; + static USOptions options = { .debug = 0, .influence_radius = DEFAULT_INFLUENCE_RADIUS_M, @@ -59,6 +75,7 @@ static USOptions options = { .lon_min = -180, .lon_max = 180, .lat_min = -90, .lat_max = 90, .cutoff_lat = 60.0 }, + .user_threads = 0, #ifdef HAVE_YAC .yac_method = -1, .yac_3d = 0, @@ -70,7 +87,6 @@ static void update_display(void); static void animation_tick(void); static void update_dim_info_current(void); static void update_dim_label(void); -static int format_time_from_units(char *out, size_t outlen, double value, const char *units); static void on_mouse_click(int px, int py); /* Callbacks */ @@ -402,10 +418,10 @@ static void on_mouse_click(int px, int py) { /* Build title */ if (current_var->units[0]) { - snprintf(ts_data.title, sizeof(ts_data.title), "%s (%s) at %.2f, %.2f", + snprintf(ts_data.title, sizeof(ts_data.title), "%.200s (%.200s) at %.2f, %.2f", current_var->name, current_var->units, lon, lat); } else { - snprintf(ts_data.title, sizeof(ts_data.title), "%s at %.2f, %.2f", + snprintf(ts_data.title, sizeof(ts_data.title), "%.200s at %.2f, %.2f", current_var->name, lon, lat); } @@ -418,7 +434,7 @@ static void on_mouse_click(int px, int py) { for (int i = 0; i < n_current_dims; i++) { if (strcmp(current_dim_info[i].name, time_dim_name) == 0) { if (current_dim_info[i].units[0]) { - strncpy(ts_data.x_label, current_dim_info[i].units, sizeof(ts_data.x_label) - 1); + snprintf(ts_data.x_label, sizeof(ts_data.x_label), "%s", current_dim_info[i].units); } break; } @@ -429,10 +445,10 @@ static void on_mouse_click(int px, int py) { } if (current_var->units[0]) { - snprintf(ts_data.y_label, sizeof(ts_data.y_label), "%s (%s)", + snprintf(ts_data.y_label, sizeof(ts_data.y_label), "%.120s (%.120s)", current_var->name, current_var->units); } else { - strncpy(ts_data.y_label, current_var->name, sizeof(ts_data.y_label) - 1); + snprintf(ts_data.y_label, sizeof(ts_data.y_label), "%s", current_var->name); } printf("Time series: %zu points (%zu valid)\n", n_out, ts_data.n_valid); @@ -523,13 +539,13 @@ static void on_range_button(void) { static void on_render_mode_toggle(void) { if (!view) return; - + /* In polygon-only mode, don't allow switching */ if (options.polygon_only) { printf("Polygon-only mode: cannot switch to interpolate mode\n"); return; } - + int result = view_toggle_render_mode(view); if (result >= 0) { const char *mode_name = (result == RENDER_MODE_POLYGON) ? "Polygon" : "Interp"; @@ -615,105 +631,6 @@ static void update_dim_info_current(void) { } } -static int parse_time_units(const char *units, double *unit_seconds, - int *y, int *mo, int *d, int *h, int *mi, double *sec) { - if (!units || !unit_seconds || !y || !mo || !d || !h || !mi || !sec) return 0; - - const char *since = strstr(units, "since"); - if (!since) return 0; - - char unit_buf[32] = {0}; - if (sscanf(units, "%31s", unit_buf) != 1) return 0; - - /* Normalize unit token to lower case */ - for (char *p = unit_buf; *p; ++p) { - if (*p >= 'A' && *p <= 'Z') *p = (char)(*p - 'A' + 'a'); - } - - if (strcmp(unit_buf, "seconds") == 0 || strcmp(unit_buf, "second") == 0 || - strcmp(unit_buf, "secs") == 0 || strcmp(unit_buf, "sec") == 0 || strcmp(unit_buf, "s") == 0) { - *unit_seconds = 1.0; - } else if (strcmp(unit_buf, "minutes") == 0 || strcmp(unit_buf, "minute") == 0 || - strcmp(unit_buf, "mins") == 0 || strcmp(unit_buf, "min") == 0) { - *unit_seconds = 60.0; - } else if (strcmp(unit_buf, "hours") == 0 || strcmp(unit_buf, "hour") == 0 || - strcmp(unit_buf, "hrs") == 0 || strcmp(unit_buf, "hr") == 0) { - *unit_seconds = 3600.0; - } else if (strcmp(unit_buf, "days") == 0 || strcmp(unit_buf, "day") == 0) { - *unit_seconds = 86400.0; - } else { - return 0; - } - - /* Parse origin date/time after "since" */ - const char *p = since + 5; - while (*p == ' ') p++; - int n = sscanf(p, "%d-%d-%d %d:%d:%lf", y, mo, d, h, mi, sec); - if (n < 3) return 0; - if (n == 3) { *h = 0; *mi = 0; *sec = 0.0; } - return 1; -} - -static int64_t days_from_civil(int y, unsigned m, unsigned d) { - y -= m <= 2; - const int era = (y >= 0 ? y : y - 399) / 400; - const unsigned yoe = (unsigned)(y - era * 400); - const unsigned doy = (153 * (m + (m > 2 ? -3 : 9)) + 2) / 5 + d - 1; - const unsigned doe = yoe * 365 + yoe / 4 - yoe / 100 + doy; - return (int64_t)(era * 146097 + (int)doe - 719468); -} - -static void civil_from_days(int64_t z, int *y, unsigned *m, unsigned *d) { - z += 719468; - const int era = (z >= 0 ? z : z - 146096) / 146097; - const unsigned doe = (unsigned)(z - era * 146097); - const unsigned yoe = (doe - doe / 1460 + doe / 36524 - doe / 146096) / 365; - int y_tmp = (int)(yoe) + era * 400; - const unsigned doy = doe - (365 * yoe + yoe / 4 - yoe / 100); - const unsigned mp = (5 * doy + 2) / 153; - const unsigned d_tmp = doy - (153 * mp + 2) / 5 + 1; - const unsigned m_tmp = mp + (mp < 10 ? 3 : -9); - y_tmp += (m_tmp <= 2); - - *y = y_tmp; - *m = m_tmp; - *d = d_tmp; -} - -static int format_time_from_units(char *out, size_t outlen, double value, const char *units) { - double unit_seconds = 0.0; - int y, mo, d, h, mi; - double sec; - if (!parse_time_units(units, &unit_seconds, &y, &mo, &d, &h, &mi, &sec)) { - return 0; - } - - int64_t days = days_from_civil(y, (unsigned)mo, (unsigned)d); - double total_sec = (double)days * 86400.0 + (double)h * 3600.0 + (double)mi * 60.0 + sec; - total_sec += value * unit_seconds; - - int64_t out_days = (int64_t)(total_sec / 86400.0); - double rem = total_sec - (double)out_days * 86400.0; - if (rem < 0) { - rem += 86400.0; - out_days -= 1; - } - - int out_y; - unsigned out_m, out_d; - civil_from_days(out_days, &out_y, &out_m, &out_d); - - int out_h = (int)(rem / 3600.0); - rem -= out_h * 3600.0; - int out_mi = (int)(rem / 60.0); - rem -= out_mi * 60.0; - int out_s = (int)(rem + 0.5); - - snprintf(out, outlen, "%04d-%02u-%02u %02d:%02d:%02d", - out_y, out_m, out_d, out_h, out_mi, out_s); - return 1; -} - static const USDimInfo *find_dim_info_for_dim(const char *dim_name) { if (!current_dim_info || !dim_name) return NULL; for (int i = 0; i < n_current_dims; i++) { @@ -874,6 +791,8 @@ static void print_usage(const char *prog) { fprintf(stderr, " --yac-3d Fractional fill-value masking for 3D variables\n"); #endif fprintf(stderr, " --light Use light theme (default: dark)\n"); + fprintf(stderr, " -t, --threads Max threads\n"); + fprintf(stderr, " (default: OMP_NUM_THREADS or 4)\n"); fprintf(stderr, " -h, --help Show this help\n"); fprintf(stderr, "\nExamples:\n"); fprintf(stderr, " %s data.nc # Single file\n", prog); @@ -882,37 +801,49 @@ static void print_usage(const char *prog) { fprintf(stderr, " %s data.1960.nc data.1961.nc -m mesh # Multi-file explicit\n", prog); } -int main(int argc, char *argv[]) { - const char *mesh_filename = NULL; - int n_data_files = 0; - const char **data_filenames = NULL; +/* + * Parse command-line options into the global `options` struct. + * Returns 0 on success, >0 for clean exit (--help), <0 on error. + * Sets *first_data_arg to the index of the first non-option argument. + */ +static int parse_options(int argc, char **argv, int *first_data_arg) { + /* Long-only option codes (above ASCII range to avoid conflicts) */ + enum { + OPT_BOX = 256, + OPT_POLAR, + OPT_CUTOFF, + OPT_YAC, + OPT_YAC_METHOD, + OPT_YAC_3D, + OPT_LIGHT, + }; - /* Parse command line options */ static struct option long_options[] = { {"mesh", required_argument, 0, 'm'}, {"resolution", required_argument, 0, 'r'}, {"influence", required_argument, 0, 'i'}, {"delay", required_argument, 0, 'd'}, {"polygon-only", no_argument, 0, 'p'}, - {"box", required_argument, 0, 1200}, - {"polar", required_argument, 0, 1201}, - {"cutoff", required_argument, 0, 1202}, + {"box", required_argument, 0, OPT_BOX}, + {"polar", required_argument, 0, OPT_POLAR}, + {"cutoff", required_argument, 0, OPT_CUTOFF}, #ifdef HAVE_YAC - {"yac", no_argument, 0, 1099}, - {"yac-method", required_argument, 0, 1100}, - {"yac-3d", no_argument, 0, 1101}, + {"yac", no_argument, 0, OPT_YAC}, + {"yac-method", required_argument, 0, OPT_YAC_METHOD}, + {"yac-3d", no_argument, 0, OPT_YAC_3D}, #endif - {"light", no_argument, 0, 1300}, + {"light", no_argument, 0, OPT_LIGHT}, + {"threads", required_argument, 0, 't'}, {"help", no_argument, 0, 'h'}, {0, 0, 0, 0} }; int opt; - while ((opt = getopt_long(argc, argv, "m:r:i:d:ph", long_options, NULL)) != -1) { + while ((opt = getopt_long(argc, argv, "m:r:i:d:pt:h", long_options, NULL)) != -1) { switch (opt) { case 'm': - mesh_filename = optarg; strncpy(options.mesh_file, optarg, MAX_NAME_LEN - 1); + options.mesh_file[MAX_NAME_LEN - 1] = '\0'; break; case 'r': options.target_resolution = atof(optarg); @@ -927,31 +858,31 @@ int main(int argc, char *argv[]) { options.polygon_only = 1; break; #ifdef HAVE_YAC - case 1099: + case OPT_YAC: options.yac_method = (int)YAC_METHOD_AVERAGE_ARITH; break; - case 1100: { + case OPT_YAC_METHOD: { USYacMethod ym; if (yac_method_parse(optarg, &ym) != 0) { fprintf(stderr, "Unknown YAC method: %s\n", optarg); fprintf(stderr, "Available: nnn1, nnn4dist, nnn4gauss, " "avg_arith, avg_dist, avg_bary, conserv1, conserv2\n"); - return 1; + return -1; } options.yac_method = (int)ym; break; } - case 1101: + case OPT_YAC_3D: options.yac_3d = 1; if (options.yac_method < 0) options.yac_method = (int)YAC_METHOD_AVERAGE_ARITH; break; #endif - case 1200: { + case OPT_BOX: { double w, e, s, n; if (sscanf(optarg, "%lf,%lf,%lf,%lf", &w, &e, &s, &n) != 4) { fprintf(stderr, "Invalid --box format, use W,E,S,N (e.g. -10,30,35,70)\n"); - return 1; + return -1; } options.target_config.lon_min = w; options.target_config.lon_max = e; @@ -959,38 +890,63 @@ int main(int argc, char *argv[]) { options.target_config.lat_max = n; break; } - case 1201: + case OPT_POLAR: if (strcmp(optarg, "north") == 0) { options.target_config.projection = PROJ_LAEA_NORTH; } else if (strcmp(optarg, "south") == 0) { options.target_config.projection = PROJ_LAEA_SOUTH; } else { fprintf(stderr, "Invalid --polar value: %s (use north or south)\n", optarg); - return 1; + return -1; } break; - case 1202: + case OPT_CUTOFF: options.target_config.cutoff_lat = atof(optarg); break; - case 1300: + case OPT_LIGHT: x_set_light_theme(); break; + case 't': { + int nt = atoi(optarg); + if (nt < 1) { + fprintf(stderr, "Invalid --threads value: %s (must be >= 1)\n", optarg); + return -1; + } + options.user_threads = nt; + break; + } case 'h': + print_usage(argv[0]); + return 1; default: print_usage(argv[0]); - return (opt == 'h') ? 0 : 1; + return -1; } } if (optind >= argc) { fprintf(stderr, "Error: No data file specified\n"); print_usage(argv[0]); - return 1; + return -1; + } + + *first_data_arg = optind; + return 0; +} + +int main(int argc, char *argv[]) { + int first_data_arg = 0; + int parse_rc = parse_options(argc, argv, &first_data_arg); + if (parse_rc != 0) { + return (parse_rc > 0) ? 0 : 1; } - /* Collect data file arguments */ - n_data_files = argc - optind; - data_filenames = (const char **)&argv[optind]; + int n_data_files = argc - first_data_arg; + const char **data_filenames = (const char **)&argv[first_data_arg]; + const char *mesh_filename = options.mesh_file[0] ? options.mesh_file : NULL; + + /* Apply thread settings: CLI > OMP_NUM_THREADS > default 4 */ + apply_thread_settings(options.user_threads); printf("=== ushow: Unstructured Data Viewer ===\n\n"); diff --git a/src/uterm.c b/src/uterm.c index 702f040..776f04d 100644 --- a/src/uterm.c +++ b/src/uterm.c @@ -5,9 +5,11 @@ * (netCDF/Zarr -> mesh -> regrid -> view) and renders frames as text. */ -#include "ushow.defines.h" +#include "us_types.h" #include "mesh.h" #include "regrid.h" +#include "kdtree.h" +#include "common.h" #ifdef HAVE_YAC #include "regrid_yac.h" #endif @@ -45,8 +47,6 @@ #define CP_LOWER_HALF_BLOCK 0x2584 #define CP_FULL_BLOCK 0x2588 -static int format_time_from_units(char *out, size_t outlen, double value, const char *units); - /* Global state */ static USFile *file = NULL; static USFileSet *fileset = NULL; @@ -77,6 +77,7 @@ typedef struct { char mesh_file[MAX_NAME_LEN]; char glyph_ramp[128]; USTargetConfig target_config; + int user_threads; /* CLI --threads value (0 = not set) */ #ifdef HAVE_YAC int yac_method; int yac_3d; @@ -95,6 +96,7 @@ static UTermOptions options = { .lon_min = -180, .lon_max = 180, .lat_min = -90, .lat_max = 90, .cutoff_lat = 60.0 }, + .user_threads = 0, #ifdef HAVE_YAC .yac_method = -1, .yac_3d = 0, @@ -256,6 +258,8 @@ static void print_usage(const char *prog) { fprintf(stderr, " conserv1, conserv2\n"); fprintf(stderr, " --yac-3d Fractional fill-value masking for 3D variables\n"); #endif + fprintf(stderr, " -t, --threads Max threads\n"); + fprintf(stderr, " (default: OMP_NUM_THREADS or 4)\n"); fprintf(stderr, " -h, --help Show this help\n\n"); fprintf(stderr, "Keys:\n"); @@ -394,84 +398,6 @@ static void save_frame(void) { } } -static int parse_time_units(const char *units, double *unit_seconds, - int *origin_year, int *origin_month, int *origin_day, - int *origin_hour, int *origin_minute, int *origin_second) { - if (!units || !unit_seconds || !origin_year || !origin_month || !origin_day || - !origin_hour || !origin_minute || !origin_second) { - return 0; - } - - const char *since = strstr(units, " since "); - if (!since) return 0; - - if (strncmp(units, "days", since - units) == 0) { - *unit_seconds = 86400.0; - } else if (strncmp(units, "day", since - units) == 0) { - *unit_seconds = 86400.0; - } else if (strncmp(units, "hours", since - units) == 0) { - *unit_seconds = 3600.0; - } else if (strncmp(units, "hour", since - units) == 0) { - *unit_seconds = 3600.0; - } else if (strncmp(units, "minutes", since - units) == 0) { - *unit_seconds = 60.0; - } else if (strncmp(units, "minute", since - units) == 0) { - *unit_seconds = 60.0; - } else if (strncmp(units, "seconds", since - units) == 0) { - *unit_seconds = 1.0; - } else if (strncmp(units, "second", since - units) == 0) { - *unit_seconds = 1.0; - } else { - return 0; - } - - const char *origin = since + strlen(" since "); - int y = 0, mo = 0, d = 0, h = 0, mi = 0, s = 0; - if (sscanf(origin, "%d-%d-%d %d:%d:%d", &y, &mo, &d, &h, &mi, &s) < 3) { - return 0; - } - - *origin_year = y; - *origin_month = mo; - *origin_day = d; - *origin_hour = h; - *origin_minute = mi; - *origin_second = s; - return 1; -} - -static int format_time_from_units(char *out, size_t outlen, double value, const char *units) { - if (!out || outlen == 0 || !units) return 0; - - double unit_seconds = 0.0; - int y = 0, mo = 0, d = 0, h = 0, mi = 0, s = 0; - if (!parse_time_units(units, &unit_seconds, &y, &mo, &d, &h, &mi, &s)) { - return 0; - } - - struct tm origin_tm = {0}; - origin_tm.tm_year = y - 1900; - origin_tm.tm_mon = mo - 1; - origin_tm.tm_mday = d; - origin_tm.tm_hour = h; - origin_tm.tm_min = mi; - origin_tm.tm_sec = s; - - time_t origin_time = timegm(&origin_tm); - if (origin_time == (time_t)-1) return 0; - - double total_seconds = value * unit_seconds; - time_t target_time = origin_time + (time_t)llround(total_seconds); - struct tm result_tm; - if (!gmtime_r(&target_time, &result_tm)) return 0; - - if (result_tm.tm_hour == 0 && result_tm.tm_min == 0 && result_tm.tm_sec == 0) { - return strftime(out, outlen, "%Y-%m-%d", &result_tm) > 0; - } - - return strftime(out, outlen, "%Y-%m-%d %H:%M:%S", &result_tm) > 0; -} - static void render_frame(int show_help, int animating) { if (!view || !current_var) return; @@ -519,7 +445,7 @@ static void render_frame(int show_help, int animating) { strncpy(time_stamp, formatted, sizeof(time_stamp) - 1); time_stamp[sizeof(time_stamp) - 1] = '\0'; } else if (current_dim_info[i].units[0]) { - snprintf(time_stamp, sizeof(time_stamp), " %.6g %s", + snprintf(time_stamp, sizeof(time_stamp), " %.6g %.44s", v, current_dim_info[i].units); } else { snprintf(time_stamp, sizeof(time_stamp), " %.6g", v); @@ -930,29 +856,44 @@ static void cleanup_all(void) { } static int parse_options(int argc, char **argv, int *first_data_arg) { + /* Long-only option codes (above ASCII range to avoid conflicts) */ + enum { + OPT_CHARS = 256, + OPT_COLOR, + OPT_NO_COLOR, + OPT_RENDER, + OPT_BOX, + OPT_POLAR, + OPT_CUTOFF, + OPT_YAC, + OPT_YAC_METHOD, + OPT_YAC_3D, + }; + static struct option long_options[] = { {"mesh", required_argument, 0, 'm'}, {"resolution", required_argument, 0, 'r'}, {"influence", required_argument, 0, 'i'}, {"delay", required_argument, 0, 'd'}, - {"chars", required_argument, 0, 1000}, - {"render", required_argument, 0, 1003}, - {"color", no_argument, 0, 1001}, - {"no-color", no_argument, 0, 1002}, - {"box", required_argument, 0, 1200}, - {"polar", required_argument, 0, 1201}, - {"cutoff", required_argument, 0, 1202}, + {"chars", required_argument, 0, OPT_CHARS}, + {"render", required_argument, 0, OPT_RENDER}, + {"color", no_argument, 0, OPT_COLOR}, + {"no-color", no_argument, 0, OPT_NO_COLOR}, + {"box", required_argument, 0, OPT_BOX}, + {"polar", required_argument, 0, OPT_POLAR}, + {"cutoff", required_argument, 0, OPT_CUTOFF}, #ifdef HAVE_YAC - {"yac", no_argument, 0, 1099}, - {"yac-method", required_argument, 0, 1100}, - {"yac-3d", no_argument, 0, 1101}, + {"yac", no_argument, 0, OPT_YAC}, + {"yac-method", required_argument, 0, OPT_YAC_METHOD}, + {"yac-3d", no_argument, 0, OPT_YAC_3D}, #endif + {"threads", required_argument, 0, 't'}, {"help", no_argument, 0, 'h'}, {0, 0, 0, 0} }; int opt; - while ((opt = getopt_long(argc, argv, "m:r:i:d:h", long_options, NULL)) != -1) { + while ((opt = getopt_long(argc, argv, "m:r:i:d:t:h", long_options, NULL)) != -1) { switch (opt) { case 'm': strncpy(options.mesh_file, optarg, MAX_NAME_LEN - 1); @@ -971,11 +912,11 @@ static int parse_options(int argc, char **argv, int *first_data_arg) { case 'h': print_usage(argv[0]); return 1; - case 1000: + case OPT_CHARS: strncpy(options.glyph_ramp, optarg, sizeof(options.glyph_ramp) - 1); options.glyph_ramp[sizeof(options.glyph_ramp) - 1] = '\0'; break; - case 1003: { + case OPT_RENDER: { int mode = TERM_RENDER_ASCII; if (term_parse_render_mode(optarg, &mode) != 0) { fprintf(stderr, "Invalid render mode: %s (use ascii|half|braille)\n", optarg); @@ -984,17 +925,17 @@ static int parse_options(int argc, char **argv, int *first_data_arg) { options.render_mode = mode; break; } - case 1001: + case OPT_COLOR: options.color_mode = 1; break; - case 1002: + case OPT_NO_COLOR: options.color_mode = 0; break; #ifdef HAVE_YAC - case 1099: + case OPT_YAC: options.yac_method = (int)YAC_METHOD_AVERAGE_ARITH; break; - case 1100: { + case OPT_YAC_METHOD: { USYacMethod ym; if (yac_method_parse(optarg, &ym) != 0) { fprintf(stderr, "Unknown YAC method: %s\n", optarg); @@ -1005,13 +946,13 @@ static int parse_options(int argc, char **argv, int *first_data_arg) { options.yac_method = (int)ym; break; } - case 1101: + case OPT_YAC_3D: options.yac_3d = 1; if (options.yac_method < 0) options.yac_method = (int)YAC_METHOD_AVERAGE_ARITH; break; #endif - case 1200: { + case OPT_BOX: { double w, e, s, n; if (sscanf(optarg, "%lf,%lf,%lf,%lf", &w, &e, &s, &n) != 4) { fprintf(stderr, "Invalid --box format, use W,E,S,N (e.g. -10,30,35,70)\n"); @@ -1023,7 +964,7 @@ static int parse_options(int argc, char **argv, int *first_data_arg) { options.target_config.lat_max = n; break; } - case 1201: + case OPT_POLAR: if (strcmp(optarg, "north") == 0) { options.target_config.projection = PROJ_LAEA_NORTH; } else if (strcmp(optarg, "south") == 0) { @@ -1033,9 +974,18 @@ static int parse_options(int argc, char **argv, int *first_data_arg) { return -1; } break; - case 1202: + case OPT_CUTOFF: options.target_config.cutoff_lat = atof(optarg); break; + case 't': { + int nt = atoi(optarg); + if (nt < 1) { + fprintf(stderr, "Invalid --threads value: %s (must be >= 1)\n", optarg); + return -1; + } + options.user_threads = nt; + break; + } default: print_usage(argv[0]); return -1; @@ -1048,6 +998,9 @@ static int parse_options(int argc, char **argv, int *first_data_arg) { return -1; } + /* Apply thread settings: CLI > OMP_NUM_THREADS > default 4 */ + apply_thread_settings(options.user_threads); + *first_data_arg = optind; return 0; } diff --git a/src/view.h b/src/view.h index 9f49ac0..e6960cd 100644 --- a/src/view.h +++ b/src/view.h @@ -5,7 +5,7 @@ #ifndef VIEW_H #define VIEW_H -#include "ushow.defines.h" +#include "us_types.h" /* * Create view for displaying data. diff --git a/tests/Makefile b/tests/Makefile index c83fb6f..019f7dd 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -14,7 +14,7 @@ NETCDF_LIBDIR := $(shell $(NC_CONFIG) --prefix 2>/dev/null)/lib NETCDF_RPATH := -Wl,-rpath,$(NETCDF_LIBDIR) CFLAGS += $(NETCDF_CFLAGS) -LIBS = $(NETCDF_LIBS) $(NETCDF_RPATH) -lm +LIBS = $(NETCDF_LIBS) $(NETCDF_RPATH) -lm -lpthread # Zarr support (optional) - build with: make WITH_ZARR=1 ifdef WITH_ZARR diff --git a/tests/test_colormaps.c b/tests/test_colormaps.c index e14ef1c..e11ceff 100644 --- a/tests/test_colormaps.c +++ b/tests/test_colormaps.c @@ -3,7 +3,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/colormaps.h" #include #include diff --git a/tests/test_file_grib.c b/tests/test_file_grib.c index b2ce541..2028e9f 100644 --- a/tests/test_file_grib.c +++ b/tests/test_file_grib.c @@ -5,7 +5,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #ifdef HAVE_GRIB diff --git a/tests/test_file_mitgcm.c b/tests/test_file_mitgcm.c index fb01205..09f0fbd 100644 --- a/tests/test_file_mitgcm.c +++ b/tests/test_file_mitgcm.c @@ -6,7 +6,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/file_mitgcm.h" #include "../src/mesh.h" #include diff --git a/tests/test_file_netcdf.c b/tests/test_file_netcdf.c index 26b409f..78bb9f7 100644 --- a/tests/test_file_netcdf.c +++ b/tests/test_file_netcdf.c @@ -4,7 +4,7 @@ #include "test_framework.h" #include "test_utils.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/file_netcdf.h" #include "../src/mesh.h" #include diff --git a/tests/test_file_zarr.c b/tests/test_file_zarr.c index 53dbd98..c369cb8 100644 --- a/tests/test_file_zarr.c +++ b/tests/test_file_zarr.c @@ -5,7 +5,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #ifdef HAVE_ZARR diff --git a/tests/test_integration.c b/tests/test_integration.c index d7a70d9..b5b31bd 100644 --- a/tests/test_integration.c +++ b/tests/test_integration.c @@ -6,7 +6,7 @@ #include "test_framework.h" #include "test_utils.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/regrid.h" #include "../src/colormaps.h" diff --git a/tests/test_latband.c b/tests/test_latband.c index c454a51..14c0904 100644 --- a/tests/test_latband.c +++ b/tests/test_latband.c @@ -9,7 +9,7 @@ #include "test_framework.h" #define _USE_MATH_DEFINES #include -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/healpix.h" #include diff --git a/tests/test_mesh.c b/tests/test_mesh.c index 38444ac..29e1016 100644 --- a/tests/test_mesh.c +++ b/tests/test_mesh.c @@ -6,7 +6,7 @@ #include "test_utils.h" #define _USE_MATH_DEFINES #include -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include diff --git a/tests/test_projection.c b/tests/test_projection.c index eb206fb..e7a5d84 100644 --- a/tests/test_projection.c +++ b/tests/test_projection.c @@ -5,7 +5,7 @@ #include "test_framework.h" #define _USE_MATH_DEFINES #include -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/regrid.h" #include "../src/projection.h" diff --git a/tests/test_regional_range.c b/tests/test_regional_range.c index 41d26c3..ee6e77a 100644 --- a/tests/test_regional_range.c +++ b/tests/test_regional_range.c @@ -8,7 +8,7 @@ #include "test_framework.h" #include "test_utils.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/regrid.h" #include "../src/file_netcdf.h" diff --git a/tests/test_regrid.c b/tests/test_regrid.c index a47c219..eff5d92 100644 --- a/tests/test_regrid.c +++ b/tests/test_regrid.c @@ -5,7 +5,7 @@ #include "test_framework.h" #define _USE_MATH_DEFINES #include -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/regrid.h" #include diff --git a/tests/test_timeseries.c b/tests/test_timeseries.c index 57858df..6265973 100644 --- a/tests/test_timeseries.c +++ b/tests/test_timeseries.c @@ -7,7 +7,7 @@ #include "test_framework.h" #include "test_utils.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/file_netcdf.h" #include "../src/mesh.h" #include diff --git a/tests/test_yac_3d.c b/tests/test_yac_3d.c index a76b024..11924cd 100644 --- a/tests/test_yac_3d.c +++ b/tests/test_yac_3d.c @@ -8,7 +8,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/regrid_yac.h" #include diff --git a/tests/test_yac_click.c b/tests/test_yac_click.c index 45b9aac..4c0efac 100644 --- a/tests/test_yac_click.c +++ b/tests/test_yac_click.c @@ -10,7 +10,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/kdtree.h" #include "../src/regrid_yac.h" diff --git a/tests/test_yac_curve2d.c b/tests/test_yac_curve2d.c index f15072d..5eb9412 100644 --- a/tests/test_yac_curve2d.c +++ b/tests/test_yac_curve2d.c @@ -9,7 +9,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/regrid_yac.h" #include diff --git a/tests/test_yac_regional.c b/tests/test_yac_regional.c index e921642..5e449eb 100644 --- a/tests/test_yac_regional.c +++ b/tests/test_yac_regional.c @@ -11,7 +11,7 @@ */ #include "test_framework.h" -#include "../src/ushow.defines.h" +#include "../src/us_types.h" #include "../src/mesh.h" #include "../src/projection.h" #include "../src/regrid_yac.h"