diff --git a/include/grass/defs/gprojects.h b/include/grass/defs/gprojects.h index 1270ed66c5a..2f21ba91230 100644 --- a/include/grass/defs/gprojects.h +++ b/include/grass/defs/gprojects.h @@ -9,6 +9,8 @@ int GPJ_transform(const struct pj_info *, const struct pj_info *, int GPJ_transform_array(const struct pj_info *, const struct pj_info *, const struct pj_info *, int, double *, double *, double *, int); +void GPJ_clone_transform(const struct pj_info *, struct gpj_transform_clone *); +void GPJ_free_transform_clone(struct gpj_transform_clone *); /* old API, to be removed */ int pj_do_proj(double *, double *, const struct pj_info *, diff --git a/include/grass/gprojects.h b/include/grass/gprojects.h index 803eecf8477..8fc71b8f2a9 100644 --- a/include/grass/gprojects.h +++ b/include/grass/gprojects.h @@ -48,6 +48,15 @@ struct pj_info { char *wkt; }; +/* Per-thread clone of a transform, filled by GPJ_clone_transform() and + * released by GPJ_free_transform_clone(). Bundles the cloned transform with + * the private PROJ context it was cloned into, so ownership is a single unit. + */ +struct gpj_transform_clone { + struct pj_info info; + PJ_CONTEXT *ctx; +}; + struct gpj_datum { char *name, *longname, *ellps; double dx, dy, dz; diff --git a/lib/proj/do_proj.c b/lib/proj/do_proj.c index 22c635310cc..e3773007ecf 100644 --- a/lib/proj/do_proj.c +++ b/lib/proj/do_proj.c @@ -1414,3 +1414,46 @@ int pj_do_transform(int count, double *x, double *y, double *h, } return ok; } + +/*! + * \brief Clone a transform into a fresh per-thread PROJ context + * + * PROJ transformation objects are not safe for concurrent use, so each thread + * needs its own. This fills \p clone with a copy of \p src whose PJ is cloned + * into a new private context. Release it with GPJ_free_transform_clone(). + * + * Safe to call concurrently from multiple threads with the same \p src, + * provided \p src is not modified during the calls. Each call clones into its + * own new context and touches no shared mutable state. + * + * \param src source transform (as set up by GPJ_init_transform()) + * \param[out] clone receives the per-thread clone (info plus private context) + */ +void GPJ_clone_transform(const struct pj_info *src, + struct gpj_transform_clone *clone) +{ + clone->ctx = proj_context_create(); + /* A failed context leaves the thread with no usable transform, so this + * aborts the run. */ + if (clone->ctx == NULL) + G_fatal_error(_("proj_context_create() failed for a per-thread " + "transform clone")); + clone->info = *src; + clone->info.pj = proj_clone(clone->ctx, src->pj); + if (clone->info.pj == NULL) + G_fatal_error(_("proj_clone() failed for a per-thread transform " + "clone")); +} + +/*! + * \brief Free a per-thread transform clone and its context + * + * \param clone clone filled by GPJ_clone_transform(); its cloned PJ is set to + * NULL after release + */ +void GPJ_free_transform_clone(struct gpj_transform_clone *clone) +{ + proj_destroy(clone->info.pj); + proj_context_destroy(clone->ctx); + clone->info.pj = NULL; +} diff --git a/raster/r.proj/benchmark/benchmark_r_proj.py b/raster/r.proj/benchmark/benchmark_r_proj.py new file mode 100644 index 00000000000..88cf48468ae --- /dev/null +++ b/raster/r.proj/benchmark/benchmark_r_proj.py @@ -0,0 +1,128 @@ +"""This is a benchmark script for r.proj thread scaling with grass.benchmark. + +This script sweeps through a raster size at a fixed memory setting, and then +memory at a fixed size. It then plots time, speedup, and efficiency. Creates +its own source raster and projects, so it runs standalone with +grass --exec python benchmark_r_proj.py.""" + +import os +import tempfile + +from grass.exceptions import CalledModuleError, GrassError +from grass.pygrass.modules import Module +import grass.script as gs +import grass.benchmark as bm + +# Baselines held fixed while one dimension is swept. +BASE_MAPSIZE = 50e6 # cells +BASE_MEMORY = 300 # MB +MAPSIZES = [10e6, 50e6, 100e6] +MEMORIES = [50, 100, 300, 1000] +METRICS = ["time", "speedup", "efficiency"] +MAX_NPROCS = 8 +REPEAT = 3 + +SRC_PROJECT = "src4326" +DST_PROJECT = "dst3857" +INPUT = "benchmark_r_proj_reference" +OUTPUT = "benchmark_r_proj" + + +def main(): + gisdbase = tempfile.mkdtemp(prefix="bench_r_proj_") + gs.create_project(os.path.join(gisdbase, SRC_PROJECT), epsg="4326") + gs.create_project(os.path.join(gisdbase, DST_PROJECT), epsg="3857") + + # Sweep raster size at the baseline memory. + results = [] + for mapsize in MAPSIZES: + benchmark( + gisdbase, + size=int(mapsize**0.5), + memory=BASE_MEMORY, + label=f"r.proj_{int(mapsize / 1e6)}M", + results=results, + ) + plot(results, "rastersize") + + # Sweep memory at the baseline raster size. + results = [] + for memory in MEMORIES: + benchmark( + gisdbase, + size=int(BASE_MAPSIZE**0.5), + memory=memory, + label=f"r.proj_memory_{memory}MB", + results=results, + ) + plot(results, "memory") + + +def benchmark(gisdbase, size, memory, label, results): + generate_input(gisdbase, size) + with gs.setup.init( + os.path.join(gisdbase, DST_PROJECT), env=os.environ.copy() + ) as session: + env = session.env + # Output region from r.proj's own suggested bounds for this input. + text = gs.read_command( + "r.proj", + project=SRC_PROJECT, + mapset="PERMANENT", + dbase=gisdbase, + input=INPUT, + method="nearest", + flags="g", + env=env, + ) + region = dict(token.split("=") for token in text.split()) + gs.run_command("g.region", env=env, **region) + + module = Module( + "r.proj", + project=SRC_PROJECT, + mapset="PERMANENT", + dbase=gisdbase, + input=INPUT, + output=OUTPUT, + method="nearest", + memory=memory, + env_=env, + run_=False, + overwrite=True, + ) + results.append( + bm.benchmark_nprocs( + module, label=label, max_nprocs=MAX_NPROCS, repeat=REPEAT + ) + ) + + +def generate_input(gisdbase, size): + """Generate the source raster in the EPSG:4326 project. Uses + r.surf.fractal, or r.random.surface when FFTW is unavailable.""" + with gs.setup.init( + os.path.join(gisdbase, SRC_PROJECT), env=os.environ.copy() + ) as session: + env = session.env + gs.run_command( + "g.region", n=50, s=40, w=-110, e=-90, rows=size, cols=size, env=env + ) + try: + Module("r.surf.fractal", output=INPUT, overwrite=True, env_=env) + except (CalledModuleError, GrassError): + Module("r.random.surface", output=INPUT, overwrite=True, env_=env) + + +def plot(results, sweep): + for metric in METRICS: + bm.nprocs_plot( + results, + filename=f"r_proj_{sweep}_{metric}.svg", + title=f"r.proj {sweep} {metric}", + metric=metric, + ) + + +if __name__ == "__main__": + main() diff --git a/raster/r.proj/footprint.c b/raster/r.proj/footprint.c new file mode 100644 index 00000000000..26fc83830d2 --- /dev/null +++ b/raster/r.proj/footprint.c @@ -0,0 +1,403 @@ +/* + * footprint.c - grid of input row spans for the output map. + * + * Each cell covers one output row and one column block and holds the range of + * input rows that block reaches. + */ + +#include +#include + +#include +#include + +#include "r.proj.h" + +/* The output is cut into nb column slices called blocks. A band is a + chunk of consecutive output rows processed together. A tile is one + band tall and k blocks wide and each tile is read as one strip. */ + +/* Geographic poles that land inside the output map, each stored as its output + * position and its input row. Empty when no pole lands inside. */ +struct pole_set { + int n; /* active poles, 0 to 2 */ + double ox[2], oy[2]; /* pole coordinates in the output CRS */ + double pole_row[2]; /* pole input row index */ +}; + +struct fp_cell { + /* lowest and highest input row one block of one output row needs */ + double rmin, rmax; /* rmax below rmin marks an empty cell */ +}; + +struct footprint { + int orows; /* output rows, one grid row per output row */ + int nb; /* output column blocks, at most 32 */ + int ocols; /* output columns */ + int irows; /* number of rows in the input map */ + struct fp_cell *cell; /* orows by nb cells in row major order */ +}; + +/* A block's input row range comes from projecting its first and last + column. Those two samples give the exact range only when the input + row changes monotonically across the block. When it does not, a + middle column can sit a fraction of a row outside the two samples. + The margin below and the strip reload in main.c cover that case. */ +#define FP_SAMPLING_MARGIN 1.0 + +/* When ocols does not divide evenly by nb, the extra columns are spread + out so no two blocks differ by more than one column. */ +/* Returns the first output column of block b. */ +static int block_first_col(const struct footprint *g, int b) +{ + return (int)((long long)b * g->ocols / g->nb); +} + +/* Returns the block that contains output column c. */ +static int block_of_col(const struct footprint *g, int c) +{ + return (int)(((long long)(c + 1) * g->nb - 1) / g->ocols); +} + +/* Projects the center of output cell (r, c) to an input row index. Returns 0 on + * a failed transform and leaves ri unchanged. */ +static int sample_row_index(const struct Cell_head *ohd, + const struct Cell_head *ihd, + const struct pj_info *oproj, + const struct pj_info *iproj, + const struct pj_info *tproj, const double *y_center, + int r, int c, double *ri) +{ + double xx = ohd->west + (c + 0.5) * ohd->ew_res; + double yy = y_center[r]; + + if (GPJ_transform(oproj, iproj, tproj, PJ_FWD, &xx, &yy, NULL) < 0 || + !isfinite(yy)) + return 0; + *ri = (ihd->north - yy) / ihd->ns_res; + return 1; +} + +/* Widens cell (r, b) to include any pole whose output point falls inside the + * cell rectangle. */ +static void fold_poles(const struct footprint *g, const struct Cell_head *ohd, + const struct pole_set *poles, int r, int b, + struct fp_cell *cell) +{ + if (!poles) + return; + + int c0 = block_first_col(g, b), c1 = block_first_col(g, b + 1), k; + double x_lo = ohd->west + c0 * ohd->ew_res; + double x_hi = ohd->west + c1 * ohd->ew_res; + double y_lo = ohd->north - (r + 1) * ohd->ns_res; + double y_hi = ohd->north - r * ohd->ns_res; + + for (k = 0; k < poles->n; k++) { + if (poles->ox[k] < x_lo || poles->ox[k] > x_hi || poles->oy[k] < y_lo || + poles->oy[k] > y_hi) + continue; + if (poles->pole_row[k] < cell->rmin) + cell->rmin = poles->pole_row[k]; + if (poles->pole_row[k] > cell->rmax) + cell->rmax = poles->pole_row[k]; + } +} + +/* For a lat/lon input, projects the north and south poles into the output and + * records each pole's input row, clamped to the map. A pole is the highest or + * lowest latitude, which the column samples can step over, so keeping its row + * makes sure the loaded strip reaches it. Leaves the set empty when no pole + * lands inside the output map. */ +static void build_pole_set(const struct Cell_head *ihd, + const struct pj_info *oproj, + const struct pj_info *iproj, + const struct pj_info *tproj, struct pole_set *poles) +{ + poles->n = 0; + if (ihd->proj != PROJECTION_LL) + return; + double polelat[2] = {90.0, -90.0}; + + for (int p = 0; p < 2; p++) { + double px = 0.0, py = polelat[p]; + + if (GPJ_transform(oproj, iproj, tproj, PJ_INV, &px, &py, NULL) < 0 || + !isfinite(px) || !isfinite(py)) + continue; + /* The pole sits outside the input's latitude range while its + projected position still ends up inside an output cell. The two edge + samples of that cell never reach the top of the input, so this line + forces the first or last input row into the cell's range. Without + this the strip for the cell holding the pole is too short. */ + double ri = (ihd->north - polelat[p]) / ihd->ns_res; + if (ri < 0) + ri = 0; + else if (ri > ihd->rows - 1) + ri = ihd->rows - 1; + poles->ox[poles->n] = px; + poles->oy[poles->n] = py; + poles->pole_row[poles->n] = ri; + poles->n++; + } +} + +/* Widens every non-empty cell by the sampling margin. */ +static void apply_sampling_margin(struct footprint *g) +{ + size_t n = (size_t)g->orows * g->nb, i; + + for (i = 0; i < n; i++) { + struct fp_cell *cell = &g->cell[i]; + + /* A cell stays empty when both of its boundary columns fail to + transform. */ + if (cell->rmax >= cell->rmin) { + cell->rmin -= FP_SAMPLING_MARGIN; + cell->rmax += FP_SAMPLING_MARGIN; + } + } +} + +/* Builds the footprint from block boundary samples, folds in any poles, and + * applies the sampling margin. */ +struct footprint *fp_create(const struct Cell_head *ohd, + const struct Cell_head *ihd, + const struct pj_info *oproj, + const struct pj_info *iproj, + const struct pj_info *tproj, const double *y_center) +{ + struct footprint *g = G_malloc(sizeof *g); + struct pole_set poles; + int r, b; + double *bnd; + + build_pole_set(ihd, oproj, iproj, tproj, &poles); + + g->orows = ohd->rows; + g->nb = MIN(ohd->cols, 32); + g->ocols = ohd->cols; + g->irows = ihd->rows; + g->cell = G_malloc(sizeof *g->cell * (size_t)g->orows * g->nb); + /* nb blocks need nb plus 1 boundary columns. */ + bnd = G_malloc(sizeof *bnd * ((size_t)g->nb + 1)); + + for (r = 0; r < g->orows; r++) { + /* Sample the number of blocks (g->nb) plus one boundary columns for + * this row. The last boundary is the final column. */ + int k; + + for (k = 0; k <= g->nb; k++) { + int c = k < g->nb ? block_first_col(g, k) : g->ocols - 1; + + if (!sample_row_index(ohd, ihd, oproj, iproj, tproj, y_center, r, c, + &bnd[k])) + bnd[k] = DBL_MAX; /* a failed sample is left out of the range */ + } + for (b = 0; b < g->nb; b++) { + struct fp_cell *cell = &g->cell[(size_t)r * g->nb + b]; + double lo = MIN(bnd[b], bnd[b + 1]); + double hi = MAX(bnd[b], bnd[b + 1]); + + cell->rmin = DBL_MAX; + cell->rmax = -DBL_MAX; + if (bnd[b] != DBL_MAX && bnd[b + 1] != DBL_MAX) { + cell->rmin = lo; + cell->rmax = hi; + } + else if (bnd[b] != DBL_MAX) { + cell->rmin = cell->rmax = bnd[b]; + } + else if (bnd[b + 1] != DBL_MAX) { + cell->rmin = cell->rmax = bnd[b + 1]; + } + fold_poles(g, ohd, &poles, r, b, cell); + } + } + G_free(bnd); + apply_sampling_margin(g); + return g; +} + +/* Fills *imin and *imax with the input rows needed for the given output + rows and columns. end_row is one past the last row and end_col is one + past the last column, same as the loop bounds, so the last band ends + with end_row equal to orows. Looks at every grid cell in that rectangle, + skips the empty ones, takes the min and max, then adds two pad rows for + the bicubic and lanczos method reads. */ +void fp_span(const struct footprint *g, int first_row, int end_row, + int first_col, int end_col, int *imin, int *imax) +{ + double rmin = DBL_MAX, rmax = -DBL_MAX; + int b_lo = block_of_col(g, first_col), b_hi = block_of_col(g, end_col - 1); + int r, b; + + if (end_row > g->orows) + G_fatal_error(_("Footprint grid has %d rows but output row %d was " + "requested"), + g->orows, end_row - 1); + + for (r = first_row; r < end_row; r++) + for (b = b_lo; b <= b_hi; b++) { + const struct fp_cell *cell = &g->cell[(size_t)r * g->nb + b]; + + if (cell->rmax < cell->rmin) + continue; /* empty cell */ + if (cell->rmin < rmin) + rmin = cell->rmin; + if (cell->rmax > rmax) + rmax = cell->rmax; + } + + if (rmax < rmin) { + /* An empty span reads as zero rows because the caller computes imax + minus imin plus one. */ + *imin = 0; + *imax = -1; + return; + } + /* Two pad rows each way keep the bicubic and lanczos neighbor reads inside + the strip. */ + int lo = (int)floor(rmin) - 2; + int hi = (int)floor(rmax) + 2; + + /* The pad can step past the first or last input row near the map edges. */ + if (lo < 0) + lo = 0; + if (hi > g->irows - 1) + hi = g->irows - 1; + *imin = lo; + *imax = hi; +} + +/* Block count for main.c's tile loop. */ +int fp_num_blocks(const struct footprint *g) +{ + return g->nb; +} + +/* First output column of block b. b equal to nb gives the output width. */ +int fp_block_start(const struct footprint *g, int b) +{ + return block_first_col(g, b); +} + +/* Largest number of input rows any k block wide tile of this band needs. */ +static int tallest_tile_rows(const struct footprint *g, int first_row, + int end_row, int k) +{ + int tallest = 0, tile_start; + + for (tile_start = 0; tile_start < g->nb; tile_start += k) { + int tile_end = MIN(tile_start + k, g->nb); + int imin, imax, rows; + + fp_span(g, first_row, end_row, block_first_col(g, tile_start), + block_first_col(g, tile_end), &imin, &imax); + rows = imax - imin + 1; + if (rows > tallest) + tallest = rows; + } + return tallest; +} + +/* Finds the widest tile, counted in whole blocks, whose input strip plus + output buffers fit under cap_bytes. Returns zero when even a single block + tile is too big. first_row and end_row are the band's rows, with end_row + one past the last. cell_size is the bytes per cell and in_cols is the + input map width. out_mult is how many output band buffers exist at once. + It is two when the previous band's write overlaps the next band's + compute, otherwise one. */ +static int tile_blocks_for_band(const struct footprint *g, int first_row, + int end_row, size_t cap_bytes, int out_mult, + int cell_size, int in_cols) +{ + size_t out_bytes = (size_t)(end_row - first_row) * g->ocols * cell_size; + int k; + + if (out_mult * out_bytes > cap_bytes) + return 0; + for (k = g->nb; k >= 1; k--) { + int worst = tallest_tile_rows(g, first_row, end_row, k); + size_t strip_bytes = + worst > 0 ? (size_t)worst * in_cols * cell_size : 0; + + if (strip_bytes + out_mult * out_bytes <= cap_bytes) + return k; + } + return 0; +} + +/* Grows the band height by doubling, preferring full-width bands and tiling + * only when even one full-width row busts the cap, and takes the last fitting + * height with its widest tile. Reports the finest tile strip the fallback + * message needs and returns zero when even one tiled row busts. */ +int fp_band_geometry(const struct footprint *g, int first_row, size_t cap_bytes, + int out_mult, int cell_size, int in_cols, + int *tile_blocks_out, int *worst_block_rows) +{ + int remaining = g->orows - first_row; + int best_h = 0, best_k = 0, try_height; + + /* Measures the single first row as a one block tile, the smallest read the + module could ever do. If even that does not fit the cap, this number is + used to tell the user how much memory the parallel path would need + before the serial tile cache takes over. */ + *worst_block_rows = tallest_tile_rows(g, first_row, first_row + 1, 1); + + /* Growing a band can only add input rows, never remove them. So once a + height is too big, every bigger height is too big as well. That means + only the first height that fails matters, and doubling finds it in a + few tries instead of counting up one at a time. The step after + remaining / 2 goes straight to remaining instead of past it. */ + for (try_height = 1;; + try_height = try_height > remaining / 2 ? remaining : try_height * 2) { + /* band_h is the band height in rows. */ + int band_h = MIN(try_height, remaining); + int worst = tallest_tile_rows(g, first_row, first_row + band_h, g->nb); + size_t strip_bytes = + worst > 0 ? (size_t)worst * in_cols * cell_size : 0; + size_t out_bytes = (size_t)band_h * g->ocols * cell_size; + + if (strip_bytes + out_mult * out_bytes > cap_bytes) + break; + best_h = band_h; + if (band_h == remaining) + break; + } + if (best_h > 0) { + /* tile_blocks is nb here, meaning one full width tile. */ + *tile_blocks_out = g->nb; + return best_h; + } + + /* One full-width row busts the cap, so grow while the exhaustive scan finds + * any fitting whole-block tile. */ + for (try_height = 1;; + try_height = try_height > remaining / 2 ? remaining : try_height * 2) { + int band_h = MIN(try_height, remaining); + int k = tile_blocks_for_band(g, first_row, first_row + band_h, + cap_bytes, out_mult, cell_size, in_cols); + + if (k == 0) + break; + best_h = band_h; + best_k = k; + if (band_h == remaining) + break; + } + if (best_h == 0) { + *tile_blocks_out = 0; + return 0; + } + *tile_blocks_out = best_k; + return best_h; +} + +void fp_free(struct footprint *g) +{ + if (!g) + return; + G_free(g->cell); + G_free(g); +} diff --git a/raster/r.proj/interp_strip.c b/raster/r.proj/interp_strip.c new file mode 100644 index 00000000000..65c39df9c01 --- /dev/null +++ b/raster/r.proj/interp_strip.c @@ -0,0 +1,257 @@ +/* + * interp_strip.c - strip versions of the resampling methods. They read the + * input rows from a strip held in memory as floats, instead of the block + * cache. Nearest is handled by strip_nearest() in main.c. + */ + +#include +#include +#include +#include +#include "r.proj.h" + +/* Reads one value from the strip. A row outside the loaded range records the + * needed row through need_lo/need_hi and returns a null, so the caller nulls + * the cell and the band reloads and recomputes. */ +static inline FCELL strip_val(const struct strip *s, int r, int c, int *need_lo, + int *need_hi) +{ + if (r < s->imin || r > s->imax) { + FCELL null_val; + + if (r < s->imin && r < *need_lo) + *need_lo = r; + if (r > s->imax && r > *need_hi) + *need_hi = r; + Rast_set_f_null_value(&null_val, 1); + return null_val; + } + return ((const FCELL *)s->data)[(size_t)(r - s->imin) * s->cols + c]; +} + +void strip_bilinear(const struct strip *s, void *obufptr, int cell_type, + double col_idx, double row_idx, struct Cell_head *incellhd, + int *need_lo, int *need_hi) +{ + int row, col, i, j; + FCELL t, u, result; + FCELL c[2][2]; + + row = (int)floor(row_idx - 0.5); + col = (int)floor(col_idx - 0.5); + + /* Full-map bounds check runs before any strip read. A sample outside the + * input map is set to NULL and returned, so strip_val is never asked for a + * row outside the map. */ + if (row < 0 || row + 1 >= incellhd->rows || col < 0 || + col + 1 >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + for (i = 0; i < 2; i++) + for (j = 0; j < 2; j++) { + const FCELL cell = strip_val(s, row + i, col + j, need_lo, need_hi); + + if (Rast_is_f_null_value(&cell)) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + c[i][j] = cell; + } + + t = col_idx - 0.5 - col; + u = row_idx - 0.5 - row; + + result = Rast_interp_bilinear(t, u, c[0][0], c[0][1], c[1][0], c[1][1]); + + Rast_set_f_value(obufptr, result, cell_type); +} + +void strip_cubic(const struct strip *s, void *obufptr, int cell_type, + double col_idx, double row_idx, struct Cell_head *incellhd, + int *need_lo, int *need_hi) +{ + int row, col, i, j; + FCELL t, u, result; + FCELL val[4]; + FCELL c[4][4]; + + row = (int)floor(row_idx - 0.5); + col = (int)floor(col_idx - 0.5); + + /* Full-map bounds check runs before any strip read. */ + if (row - 1 < 0 || row + 2 >= incellhd->rows || col - 1 < 0 || + col + 2 >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + for (i = 0; i < 4; i++) + for (j = 0; j < 4; j++) { + const FCELL cell = + strip_val(s, row - 1 + i, col - 1 + j, need_lo, need_hi); + + if (Rast_is_f_null_value(&cell)) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + c[i][j] = cell; + } + + t = col_idx - 0.5 - col; + u = row_idx - 0.5 - row; + + for (i = 0; i < 4; i++) { + const FCELL *tmp = c[i]; + + val[i] = Rast_interp_cubic(t, tmp[0], tmp[1], tmp[2], tmp[3]); + } + + result = Rast_interp_cubic(u, val[0], val[1], val[2], val[3]); + + Rast_set_f_value(obufptr, result, cell_type); +} + +void strip_lanczos(const struct strip *s, void *obufptr, int cell_type, + double col_idx, double row_idx, struct Cell_head *incellhd, + int *need_lo, int *need_hi) +{ + int row, col, i, j, k; + double t, u; + FCELL result; + DCELL c[25]; + + row = (int)floor(row_idx); + col = (int)floor(col_idx); + + /* Full-map bounds check runs before any strip read. */ + if (row - 2 < 0 || row + 2 >= incellhd->rows || col - 2 < 0 || + col + 2 >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + k = 0; + for (i = 0; i < 5; i++) { + for (j = 0; j < 5; j++) { + const FCELL cell = + strip_val(s, row - 2 + i, col - 2 + j, need_lo, need_hi); + + if (Rast_is_f_null_value(&cell)) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + c[k++] = cell; + } + } + + t = col_idx - 0.5 - col; + u = row_idx - 0.5 - row; + + result = Rast_interp_lanczos(t, u, c); + + Rast_set_f_value(obufptr, result, cell_type); +} + +void strip_bilinear_f(const struct strip *s, void *obufptr, int cell_type, + double col_idx, double row_idx, + struct Cell_head *incellhd, int *need_lo, int *need_hi) +{ + int row, col; + FCELL cell; + + row = (int)floor(row_idx); + col = (int)floor(col_idx); + + if (row < 0 || row >= incellhd->rows || col < 0 || col >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + cell = strip_val(s, row, col, need_lo, need_hi); + /* if nearest is null, all the other interps will be null */ + if (Rast_is_f_null_value(&cell)) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + strip_bilinear(s, obufptr, cell_type, col_idx, row_idx, incellhd, need_lo, + need_hi); + /* fallback to nearest if bilinear is null */ + if (Rast_is_f_null_value(obufptr)) + Rast_set_f_value(obufptr, cell, cell_type); +} + +void strip_cubic_f(const struct strip *s, void *obufptr, int cell_type, + double col_idx, double row_idx, struct Cell_head *incellhd, + int *need_lo, int *need_hi) +{ + int row, col; + FCELL cell; + + row = (int)floor(row_idx); + col = (int)floor(col_idx); + + if (row < 0 || row >= incellhd->rows || col < 0 || col >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + cell = strip_val(s, row, col, need_lo, need_hi); + /* if nearest is null, all the other interps will be null */ + if (Rast_is_f_null_value(&cell)) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + strip_cubic(s, obufptr, cell_type, col_idx, row_idx, incellhd, need_lo, + need_hi); + /* fallback to bilinear if cubic is null */ + if (Rast_is_f_null_value(obufptr)) { + strip_bilinear(s, obufptr, cell_type, col_idx, row_idx, incellhd, + need_lo, need_hi); + /* fallback to nearest if bilinear is null */ + if (Rast_is_f_null_value(obufptr)) + Rast_set_f_value(obufptr, cell, cell_type); + } +} + +void strip_lanczos_f(const struct strip *s, void *obufptr, int cell_type, + double col_idx, double row_idx, struct Cell_head *incellhd, + int *need_lo, int *need_hi) +{ + int row, col; + FCELL cell; + + row = (int)floor(row_idx); + col = (int)floor(col_idx); + + if (row < 0 || row >= incellhd->rows || col < 0 || col >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + cell = strip_val(s, row, col, need_lo, need_hi); + /* if nearest is null, all the other interps will be null */ + if (Rast_is_f_null_value(&cell)) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + strip_lanczos(s, obufptr, cell_type, col_idx, row_idx, incellhd, need_lo, + need_hi); + /* fallback to bicubic if lanczos is null */ + if (Rast_is_f_null_value(obufptr)) { + strip_cubic(s, obufptr, cell_type, col_idx, row_idx, incellhd, need_lo, + need_hi); + /* fallback to bilinear if cubic is null */ + if (Rast_is_f_null_value(obufptr)) { + strip_bilinear(s, obufptr, cell_type, col_idx, row_idx, incellhd, + need_lo, need_hi); + /* fallback to nearest if bilinear is null */ + if (Rast_is_f_null_value(obufptr)) + Rast_set_f_value(obufptr, cell, cell_type); + } + } +} diff --git a/raster/r.proj/main.c b/raster/r.proj/main.c index 59702cdfc92..3a415175e67 100644 --- a/raster/r.proj/main.c +++ b/raster/r.proj/main.c @@ -66,6 +66,19 @@ #include #include "r.proj.h" +#ifdef _OPENMP +#include +static inline double rproj_wtime(void) +{ + return omp_get_wtime(); +} +#else +static inline double rproj_wtime(void) +{ + return 0.0; +} +#endif + /* modify this table to add new methods */ struct menu menu[] = { {p_nearest, "nearest", "nearest neighbor"}, @@ -77,9 +90,146 @@ struct menu menu[] = { {p_lanczos_f, "lanczos_f", "lanczos filter with fallback"}, {NULL, NULL, NULL}}; +/* menu ends with a NULL row, so it has one more slot than methods. */ +enum { N_METHODS = sizeof menu / sizeof *menu - 1 }; + static char *make_ipol_list(void); static char *make_ipol_desc(void); +/* Round whole numbers through 32 bit float so every read path + returns the same values. */ +static void quantize_cell_row(void *row, int cols, int cell_type) +{ + int i; + + if (cell_type == CELL_TYPE) { + CELL *p = row; + + for (i = 0; i < cols; i++) + if (!Rast_is_c_null_value(&p[i])) + Rast_set_f_value(&p[i], (FCELL)p[i], CELL_TYPE); + } + else if (cell_type == DCELL_TYPE) { + DCELL *p = row; + + for (i = 0; i < cols; i++) + if (!Rast_is_d_null_value(&p[i])) + Rast_set_f_value(&p[i], (FCELL)p[i], DCELL_TYPE); + } +} + +/* Nearest read from the in-RAM input strip covering rows imin to imax. Indices + * are full-map and addressed relative to imin. */ +static void strip_nearest(const struct strip *s, void *obufptr, int cell_type, + double col_idx, double row_idx, + struct Cell_head *incellhd, int *need_lo, + int *need_hi) +{ + int c = (int)floor(col_idx); + int r = (int)floor(row_idx); + int cell_size = Rast_cell_size(cell_type); + + /* A sample outside the input map is a legitimate NULL, like p_nearest. */ + if (r < 0 || r >= incellhd->rows || c < 0 || c >= incellhd->cols) { + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + /* A row inside the map but outside the loaded strip records the needed row + * and nulls the cell, so the band reloads and recomputes. */ + if (r < s->imin || r > s->imax) { + if (r < s->imin && r < *need_lo) + *need_lo = r; + if (r > s->imax && r > *need_hi) + *need_hi = r; + Rast_set_null_value(obufptr, 1, cell_type); + return; + } + + unsigned char *src = (unsigned char *)s->data + + (((size_t)(r - s->imin) * s->cols + c) * cell_size); + memcpy(obufptr, src, cell_size); +} + +/* Strip kernels in menu[] order, so slot i matches menu[i].method. */ +static const strip_func strip_kernels[] = { + strip_nearest, strip_bilinear, strip_cubic, strip_lanczos, + strip_bilinear_f, strip_cubic_f, strip_lanczos_f}; + +/* The two tables must stay in step, one strip kernel per menu entry. */ +_Static_assert(sizeof strip_kernels / sizeof *strip_kernels == N_METHODS, + "strip_kernels and menu must have one entry per method"); + +/* Grid that sizes band heights and column tiles. */ +static struct footprint *band_grid = NULL; + +/* Serial tile-cache path for output rows whose input footprint is too tall to + * band. Finishes the run from row obr0 with the readcell cache so it matches + * serial r.proj. */ +static void +fallback_serial_cache(int fdi, int fdo, int cell_type, int method, + const struct pj_info *oproj, const struct pj_info *iproj, + const struct pj_info *tproj, struct Cell_head *incellhd, + struct Cell_head *outcellhd, const double *y_center, + int obr0, const char *memory) +{ + /* readcell reads the input map, so it runs in the input env. */ + G_switch_env(); /* -> input */ + struct cache *ibuffer = readcell(fdi, memory); + G_switch_env(); /* -> output */ + func interpolate = menu[method].method; + void *obuffer = Rast_allocate_output_buf(cell_type); + int cell_size = Rast_cell_size(cell_type); + double local_x_start = outcellhd->west + (outcellhd->ew_res / 2); + + for (int row = obr0; row < outcellhd->rows; row++) { + G_percent(row - obr0, outcellhd->rows - obr0, 5); + for (int col = 0; col < outcellhd->cols; col++) { + void *obufptr = (unsigned char *)obuffer + (size_t)col * cell_size; + double x1 = local_x_start + col * outcellhd->ew_res; + double y1 = y_center[row]; + + if (GPJ_transform(oproj, iproj, tproj, PJ_FWD, &x1, &y1, NULL) < + 0) { + Rast_set_null_value(obufptr, 1, cell_type); + } + else { + double col_idx = (x1 - incellhd->west) / incellhd->ew_res; + double row_idx = (incellhd->north - y1) / incellhd->ns_res; + + interpolate(ibuffer, obufptr, cell_type, col_idx, row_idx, + incellhd); + } + } + Rast_put_row(fdo, obuffer, cell_type); + } + + release_cache(ibuffer); + G_free(obuffer); +} + +/* Write the previous band's output rows that were held for the overlapped + * write, then free the buffer. Returns when nothing is held. */ +static void flush_pending_band(int fdo, int cell_type, int cols, int cell_size, + void **pending, int r0, int r1) +{ + if (*pending == NULL) + return; + for (int wr = r0; wr < r1; wr++) + Rast_put_row(fdo, + (unsigned char *)*pending + + (size_t)(wr - r0) * cols * cell_size, + cell_type); + G_free(*pending); + *pending = NULL; +} + +/* Return the compute thread count, with nprocs overriding OMP_NUM_THREADS. */ +static int compute_nprocs(struct Option *nprocs) +{ + return G_set_omp_num_threads(nprocs); +} + int main(int argc, char **argv) { char *mapname, /* ptr to name of output layer */ @@ -98,14 +248,7 @@ int main(int argc, char **argv) overwrite, /* Overwrite */ curr_proj; /* output projection (see gis.h) */ - void *obuffer; /* buffer that holds one output row */ - - struct cache *ibuffer; /* buffer that holds the input map */ - func interpolate; /* interpolation routine */ - - double xcoord2, /* temporary x coordinates */ - ycoord2, /* temporary y coordinates */ - onorth, osouth, /* save original border coords */ + double onorth, osouth, /* save original border coords */ oeast, owest, inorth, isouth, ieast, iwest; char north_str[30], south_str[30], east_str[30], west_str[30]; @@ -134,6 +277,7 @@ int main(int argc, char **argv) *indbase, /* name of input database */ *interpol, /* interpolation method */ *memory, /* amount of memory for cache */ + *nprocs, /* number of compute threads */ *res, /* resolution of target map */ *format; /* output format */ @@ -191,6 +335,13 @@ int main(int argc, char **argv) memory = G_define_standard_option(G_OPT_MEMORYMB); + nprocs = G_define_standard_option(G_OPT_M_NPROCS); + nprocs->description = _( + "Number of threads for parallel computing. 0, the default, uses the " + "OpenMP default and honors OMP_NUM_THREADS if set. A value above zero " + "overrides OMP_NUM_THREADS, and a value below zero leaves that many " + "cores free"); + res = G_define_option(); res->key = "resolution"; res->type = TYPE_DOUBLE; @@ -276,7 +427,9 @@ int main(int argc, char **argv) if (!ipolname) G_fatal_error(_("<%s=%s> unknown %s"), interpol->key, interpol->answer, interpol->key); - interpolate = menu[method].method; + + /* Resolve the strip kernel once; menu[] and strip_kernels[] share order. */ + strip_func interp = strip_kernels[method]; mapname = outmap->answer ? outmap->answer : inmap->answer; if (mapname && !list->answer && !overwrite && !print_bounds->answer && @@ -420,7 +573,7 @@ int main(int argc, char **argv) if (G_verbose() > G_verbose_std()) pj_print_proj_params(&iproj, &oproj); - /* this call causes r.proj to read the entire map into memory */ + /* Read the input map's rows, columns, resolution, and bounds. */ Rast_get_cellhd(inmap->answer, setname, &incellhd); if (G_projection() == PROJECTION_XY) @@ -664,19 +817,34 @@ int main(int argc, char **argv) G_message(_("NS-res: %f"), outcellhd.ns_res); G_message(" "); - /* open and read the relevant parts of the input map and close it */ + /* Open the input map in the input env. fdi stays open across the band loop + * because each band reads only its own input strip. */ G_switch_env(); Rast_set_input_window(&incellhd); fdi = Rast_open_old(inmap->answer, setname); cell_type = Rast_get_map_type(fdi); - ibuffer = readcell(fdi, memory->answer); - Rast_close(fdi); + if (strcmp(interpol->answer, "nearest") != 0) + cell_type = FCELL_TYPE; + cell_size = Rast_cell_size(cell_type); - /* And switch back to original location */ + /* Pick the read-thread count in the input env so the mask guard checks the + * source mapset's mask. Rast_disable_omp_on_mask forces serial under a mask + * or without OpenMP. Above one thread each thread opens its own fd and fdi + * serves only the serial fallback. */ + int want_nprocs = compute_nprocs(nprocs); + int read_nprocs = Rast_disable_omp_on_mask(want_nprocs); + int *fd_read = NULL; + if (read_nprocs > 1) { + fd_read = G_malloc((size_t)read_nprocs * sizeof(int)); + for (int t = 0; t < read_nprocs; t++) + fd_read[t] = Rast_open_old(inmap->answer, setname); + } + + /* Back in the output env, set the output window, init the transform, and + * open the output map. Both fds stay open and their windows survive env + * switches, so reads and writes use the right window throughout. */ G_switch_env(); Rast_set_output_window(&outcellhd); - - /* reproject from output to input */ G_unset_window(); G_set_window(&outcellhd); tproj.def = NULL; @@ -687,73 +855,410 @@ int main(int argc, char **argv) if (GPJ_init_transform(&oproj, &iproj, &tproj) < 0) G_fatal_error(_("Unable to initialize coordinate transformation")); - if (strcmp(interpol->answer, "nearest") == 0) { + if (strcmp(interpol->answer, "nearest") == 0) fdo = Rast_open_new(mapname, cell_type); - obuffer = (CELL *)Rast_allocate_output_buf(cell_type); - } - else { + else fdo = Rast_open_fp_new(mapname); - cell_type = FCELL_TYPE; - obuffer = (FCELL *)Rast_allocate_output_buf(cell_type); - } - - cell_size = Rast_cell_size(cell_type); - - xcoord2 = outcellhd.west + (outcellhd.ew_res / 2); - ycoord2 = outcellhd.north - (outcellhd.ns_res / 2); - - G_important_message(_("Projecting...")); - for (row = 0; row < outcellhd.rows; row++) { - /* obufptr = obuffer */; - G_percent(row, outcellhd.rows - 1, 2); + /* Banding runs a serial band loop of serial strip load, parallel compute + * into a per-band buffer, and an in-order per-band write. This bounds peak + * memory by the cap rather than the whole input map. */ + double cap_mb = atof(memory->answer); + size_t cap_bytes = (size_t)(cap_mb * 1024.0 * 1024.0); + /* Output write time that overlap the compute time are counted in the + * compute time. The write time only counts the writes that run on their + * own. The fallback path is not timed. */ + double t_size = 0.0, t_fill = 0.0, t_compute = 0.0, t_write = 0.0; + int n_bands = 0; + int max_tiles = 1; /* most column tiles used by any single band */ + + /* Holds the north-to-south coordinate of the center of each output row. It + * starts at the top edge and steps down one cell per row. The step uses + * repeated subtraction to match serial r.proj exactly. */ + double *y_center = G_malloc((size_t)outcellhd.rows * sizeof(double)); + { + double yc = outcellhd.north - (outcellhd.ns_res / 2); + for (int r = 0; r < outcellhd.rows; r++) { + y_center[r] = yc; + yc -= outcellhd.ns_res; + } + } -#if 0 - /* parallelization does not always work, - * segfaults in the interpolation functions - * can happen */ -#pragma omp parallel for schedule(static) + /* Build the footprint that sizes the bands. It also projects any poles into + * the output and applies the sampling margin. */ + band_grid = + fp_create(&outcellhd, &incellhd, &oproj, &iproj, &tproj, y_center); + + G_important_message(_("Projecting (banded, per-thread PROJ context)...")); + + int used_fallback = 0; /* set when the serial tile-cache fallback runs */ + int force_tilecache = getenv("R_PROJ_FORCE_TILECACHE") != NULL; + /* For tests only. FG_SHRINK_SPAN=N loads each band N rows too short at the + * top and bottom, so the reload path runs. Zero leaves the strip full + * size. */ + const char *shrink_env = getenv("FG_SHRINK_SPAN"); + int shrink_span = shrink_env ? atoi(shrink_env) : 0; + /* Rolling input-strip window. win holds input rows [win_imin, win_imax]. + * win_imax < win_imin marks it empty and forces a full read. win_cap is + * its allocated byte size, and win is freed once at fallback_done. */ + unsigned char *win = NULL; + size_t win_cap = 0; + int win_imin = 0, win_imax = -1; + /* Turn on output double buffering when more than one compute thread runs, + * so one thread can write the previous band while the next one computes. + * out_mult then reserves two output bands so the memory budget matches the + * writer. */ + int write_overlap = want_nprocs > 1; + int out_mult = write_overlap ? 2 : 1; + /* Previous band's output buffer, written by one thread while the next band + * computes. It is NULL when nothing is pending and holds rows + * [pending_r0, pending_r1). */ + void *pending_out = NULL; + int pending_r0 = 0, pending_r1 = 0; + int obr0 = 0; + while (obr0 < outcellhd.rows) { + /* Co-size this band's height and tile width from the grid. A zero + * height means even one row busts the cap so the run finishes from obr0 + * on the serial tile-cache path. */ + double ts = rproj_wtime(); + int imin = 0, imax = -1; + int tile_blocks = 0, worst_block_rows = 0; + int band_orows = + force_tilecache + ? 0 + : fp_band_geometry(band_grid, obr0, cap_bytes, out_mult, + cell_size, incellhd.cols, &tile_blocks, + &worst_block_rows); + if (band_orows == 0) { + if (force_tilecache) { + G_warning(_("R_PROJ_FORCE_TILECACHE is set: taking the serial " + "tile-cache path for all output rows (testing " + "override).")); + } + else { + size_t out1 = (size_t)outcellhd.cols * cell_size; + size_t strip_bytes = + worst_block_rows > 0 + ? (size_t)worst_block_rows * incellhd.cols * cell_size + : 0; + int needed_mb = (int)ceil((double)(strip_bytes + out1) / + (1024.0 * 1024.0)) + + 1; + G_warning( + _("Memory cap (%.1f MB) is below what one output row " + "needs (input footprint %d rows, %.1f MB). Falling " + "back to the serial tile-cache path for output rows " + "%d-%d; this path is slower. Raise memory= to at " + "least %d MB to use the parallel path."), + cap_mb, worst_block_rows, + (double)(strip_bytes + out1) / (1024.0 * 1024.0), obr0, + outcellhd.rows - 1, needed_mb); + } + /* Flush the deferred band before the fallback writes from obr0 in + * order. */ + flush_pending_band(fdo, cell_type, outcellhd.cols, cell_size, + &pending_out, pending_r0, pending_r1); + fallback_serial_cache(fdi, fdo, cell_type, method, &oproj, &iproj, + &tproj, &incellhd, &outcellhd, y_center, obr0, + memory->answer); + used_fallback = 1; + goto fallback_done; + } + t_size += rproj_wtime() - ts; + + int obr1 = obr0 + band_orows; + n_bands++; + int nb = fp_num_blocks(band_grid); + int n_tiles = (nb + tile_blocks - 1) / tile_blocks; + if (n_tiles > max_tiles) + max_tiles = n_tiles; + + /* Per-band output buffer at full width, filled tile by tile and written + * once after all tiles. The row slots are disjoint, so compute is + * lock-free. */ + void *band_out = + G_malloc((size_t)band_orows * outcellhd.cols * cell_size); + + /* Column tiles are processed one at a time, so peak strip memory is the + * worst tile rather than the band's union. A single tile spanning every + * block is the full-width fast path. */ + for (int tb = 0; tb < nb; tb += tile_blocks) { + int te = tb + tile_blocks < nb ? tb + tile_blocks : nb; + int obc0 = fp_block_start(band_grid, tb); + int obc1 = fp_block_start(band_grid, te); + + /* Fill spans come from the grid. The strip is full input width + * because the raster API reads whole rows, so columns are not + * cropped. */ + fp_span(band_grid, obr0, obr1, obc0, obc1, &imin, &imax); + /* The test hook shortens the span so the reload path runs. */ + if (shrink_span > 0 && imax >= imin) { + imin += shrink_span; + imax -= shrink_span; + } + int reloaded = 0; /* set once the span is widened and reread */ + + reload_tile:; + int strip_rows = imax - imin + 1; + + /* Serial strip load, since a single fd makes get_row unsafe to + * share. An empty tile with strip_rows at or below zero is not + * read. Its cells project outside the input and the kernel sets + * them NULL, and the window is invalidated so the next band + * re-reads in full. */ + void *strip = NULL; + if (strip_rows > 0) { + size_t need = (size_t)strip_rows * incellhd.cols * cell_size; + size_t row_bytes = (size_t)incellhd.cols * cell_size; + /* Slide only for a serial-read single-tile band whose rows + * advance forward and still overlap the window. Anything else + * does a full read, same as before. */ + int can_slide = read_nprocs == 1 && n_tiles == 1 && + win_imax >= 0 && imin >= win_imin && + imin <= win_imax + 1; + int read_from = imin; + /* Grow, then memmove, then read the tail. G_realloc may move + * the buffer, so it must run before the memmove that + * repositions the retained overlap. Realloc preserves the old + * rows at their old offsets and the memmove shifts them to the + * new imin origin. Rows past a shrunk imax are dropped through + * win_imax below and re-read if a later band needs them. */ + if (need > win_cap) { + win = G_realloc(win, need); + win_cap = need; + } + if (can_slide && win_imax >= imin) { + memmove(win, win + (size_t)(imin - win_imin) * row_bytes, + (size_t)(win_imax - imin + 1) * row_bytes); + read_from = win_imax + 1; /* only new rows hit disk */ + } + strip = win; + double t0 = rproj_wtime(); + if (read_from <= imax) { + G_switch_env(); /* -> input */ + if (read_nprocs > 1) { +#ifdef _OPENMP + /* Parallel read of the full span only, since can_slide + * is false here. Each thread reads a contiguous + * disjoint block of rows through its own fd into its + * own strip slice, so no two threads share an fd or a + * row. */ +#pragma omp parallel num_threads(read_nprocs) + { + int t = omp_get_thread_num(); +#pragma omp for schedule(static) + for (int r = read_from; r <= imax; r++) + Rast_get_row(fd_read[t], + (unsigned char *)strip + + (size_t)(r - imin) * + incellhd.cols * cell_size, + r, cell_type); + } #endif - - for (col = 0; col < outcellhd.cols; col++) { - void *obufptr = - (void *)((const unsigned char *)obuffer + col * cell_size); - - double xcoord1 = xcoord2 + (col)*outcellhd.ew_res; - double ycoord1 = ycoord2; - - /* project coordinates in output matrix to */ - /* coordinates in input matrix */ - if (GPJ_transform(&oproj, &iproj, &tproj, PJ_FWD, &xcoord1, - &ycoord1, NULL) < 0) { - G_fatal_error(_("Error in %s"), "GPJ_transform()"); - Rast_set_null_value(obufptr, 1, cell_type); + } + else { + /* Serial read of the (possibly partial) tail. */ + for (int r = read_from; r <= imax; r++) + Rast_get_row(fdi, + (unsigned char *)strip + + (size_t)(r - imin) * + incellhd.cols * cell_size, + r, cell_type); + } + G_switch_env(); /* -> output */ + for (int r = read_from; r <= imax; r++) + quantize_cell_row((unsigned char *)strip + + (size_t)(r - imin) * + incellhd.cols * cell_size, + incellhd.cols, cell_type); + } + t_fill += rproj_wtime() - t0; + /* Record what the window now holds. A tiled band leaves win + * with only its last tile, so invalidate both fields to force + * the next band's full read. Setting both keeps validity + * independent of the && short-circuit order. */ + if (n_tiles == 1) { + win_imin = imin; + win_imax = imax; + } + else { + win_imin = 0; + win_imax = -1; + } } else { - /* convert to row/column indices of input matrix */ - - /* column index in input matrix */ - double col_idx = (xcoord1 - incellhd.west) / incellhd.ew_res; - - /* row index in input matrix */ - double row_idx = (incellhd.north - ycoord1) / incellhd.ns_res; - - /* and resample data point */ - interpolate(ibuffer, obufptr, cell_type, col_idx, row_idx, - &incellhd); + win_imin = 0; + win_imax = -1; /* empty tile, nothing resident */ } - /* obufptr = G_incr_void_ptr(obufptr, cell_size); */ + struct strip sd = {strip, imin, imax, incellhd.cols}; + /* Widest input rows a kernel wants below imin and above imax. The + * reduction leaves them at imin and imax when the strip covered + * every sample. */ + int need_lo = imin, need_hi = imax; + double t1 = rproj_wtime(); + /* One parallel region per tile. The omp for divides the band's + * output rows among this region's threads. The directives are + * separate so each thread clones its PROJ context before the row + * loop and destroys it after. */ +#pragma omp parallel + { + struct gpj_transform_clone tproj_local; + GPJ_clone_transform(&tproj, &tproj_local); + +#pragma omp single nowait + { + /* One thread writes the previous band's rows in order while + * the rest compute this band, then frees it. Only the first + * tile has a pending band, and clearing it keeps a reload + * from writing the same rows twice. */ + if (obc0 == 0 && pending_out) { + for (int wr = pending_r0; wr < pending_r1; wr++) + Rast_put_row(fdo, + (unsigned char *)pending_out + + (size_t)(wr - pending_r0) * + outcellhd.cols * cell_size, + cell_type); + G_free(pending_out); + pending_out = NULL; + } + } + +#pragma omp for private(row, col) schedule(dynamic) reduction(min : need_lo) \ + reduction(max : need_hi) + for (row = obr0; row < obr1; row++) { + void *out_row = + (unsigned char *)band_out + + (size_t)(row - obr0) * outcellhd.cols * cell_size; + double local_y = y_center[row]; + double local_x_start = + outcellhd.west + (outcellhd.ew_res / 2); + + for (col = obc0; col < obc1; col++) { + void *obufptr = + (unsigned char *)out_row + (size_t)col * cell_size; + double x1 = local_x_start + (col * outcellhd.ew_res); + double y1 = local_y; + + if (GPJ_transform(&oproj, &iproj, &tproj_local.info, + PJ_FWD, &x1, &y1, NULL) < 0) { + Rast_set_null_value(obufptr, 1, cell_type); + } + else { + double c_idx = + (x1 - incellhd.west) / incellhd.ew_res; + double r_idx = + (incellhd.north - y1) / incellhd.ns_res; + interp(&sd, obufptr, cell_type, c_idx, r_idx, + &incellhd, &need_lo, &need_hi); + } + } + } + + GPJ_free_transform_clone(&tproj_local); + } + t_compute += rproj_wtime() - t1; + /* strip aliases the persistent window buffer win, so it is not + * freed per tile. win is freed once at fallback_done. */ + + /* Reload once and recompute if a kernel needed an input row the + * strip left out. */ + if (need_lo < imin || need_hi > imax) { + if (reloaded) + G_fatal_error(_("Band strip still short after a reload, " + "input rows [%d, %d] but [%d, %d] needed"), + imin, imax, need_lo, need_hi); + size_t widened = + (size_t)(need_hi - need_lo + 1) * incellhd.cols * cell_size; + size_t out_bytes = + (size_t)band_orows * outcellhd.cols * cell_size; + /* If the widened strip no longer fits the memory cap, finish + * from obr0 on the serial tile-cache path, the same path the + * band sizing falls back to. */ + if (widened + (size_t)out_mult * out_bytes > cap_bytes) { + G_free(band_out); + flush_pending_band(fdo, cell_type, outcellhd.cols, + cell_size, &pending_out, pending_r0, + pending_r1); + fallback_serial_cache(fdi, fdo, cell_type, method, &oproj, + &iproj, &tproj, &incellhd, &outcellhd, + y_center, obr0, memory->answer); + used_fallback = 1; + goto fallback_done; + } + G_verbose_message(_("Reloading band strip for output rows " + "%d-%d, input rows [%d, %d] widen to " + "[%d, %d]"), + obr0, obr1 - 1, imin, imax, need_lo, need_hi); + imin = need_lo; + imax = need_hi; + reloaded = 1; + goto reload_tile; + } } - Rast_put_row(fdo, obuffer, cell_type); + /* Defer this band so the next band's compute region writes it through + * the omp single above. The previous pending band was written and freed + * there, so pending_out is already clear. Non-overlap bands write in + * order below. */ + if (write_overlap) { + pending_out = band_out; + pending_r0 = obr0; + pending_r1 = obr1; + } + else { + double t2 = rproj_wtime(); + for (row = obr0; row < obr1; row++) + Rast_put_row(fdo, + (unsigned char *)band_out + (size_t)(row - obr0) * + outcellhd.cols * + cell_size, + cell_type); + t_write += rproj_wtime() - t2; + G_free(band_out); + } - xcoord2 = outcellhd.west + (outcellhd.ew_res / 2); - ycoord2 -= outcellhd.ns_res; + G_percent(obr1, outcellhd.rows, 5); + obr0 = obr1; } +fallback_done: + /* Flush the last band's deferred write on normal completion, timed into + * t_write. The fallback bails already flushed, so pending_out is NULL on + * those paths. */ + { + double tw = rproj_wtime(); + flush_pending_band(fdo, cell_type, outcellhd.cols, cell_size, + &pending_out, pending_r0, pending_r1); + t_write += rproj_wtime() - tw; + } + G_free(y_center); + if (band_grid) + fp_free(band_grid); + /* Single free site for the rolling window. Normal completion and both + * fallback_done bails converge here, so one free covers every path. win is + * NULL when a bail fired before any band allocated it. */ + if (win) + G_free(win); + + if (used_fallback) + G_debug(1, "PHASE_TIMERS fallback=1 fallback_from_row=%d", obr0); + else + G_debug(1, + "PHASE_TIMERS size=%.4f fill=%.4f compute=%.4f write=%.4f " + "bands=%d tiles=%d", + t_size, t_fill, t_compute, t_write, n_bands, max_tiles); + + /* Close input map in its own env, then the output map. */ + G_switch_env(); /* -> input */ + Rast_close(fdi); + if (fd_read) { + for (int t = 0; t < read_nprocs; t++) + Rast_close(fd_read[t]); + G_free(fd_read); + } + G_switch_env(); /* -> output */ Rast_close(fdo); - release_cache(ibuffer); if (have_colors > 0) { Rast_write_colors(mapname, G_mapset(), &colr); diff --git a/raster/r.proj/r.proj.h b/raster/r.proj/r.proj.h index 935415c1e41..1f6848d71d4 100644 --- a/raster/r.proj/r.proj.h +++ b/raster/r.proj/r.proj.h @@ -27,6 +27,21 @@ struct cache { typedef void (*func)(struct cache *, void *, int, double, double, struct Cell_head *); +/* One input strip resident in memory, addressed by full-map input row. Rows + * outside imin to imax are not loaded. A kernel that needs one records it + * through the need arguments so the band reloads and recomputes. */ +struct strip { + void *data; /* input rows imin to imax at full input width */ + int imin, imax; /* loaded input row range */ + int cols; /* input columns, the row stride */ +}; + +/* Strip interpolation kernels (interp_strip.c) read from an in-memory strip. + * The last two arguments carry the widest input rows a kernel wanted below imin + * and above imax, left unchanged when the strip covered every sample. */ +typedef void (*strip_func)(const struct strip *, void *, int, double, double, + struct Cell_head *, int *, int *); + struct menu { func method; /* routine to interpolate new value */ char *name; /* method name */ @@ -35,6 +50,23 @@ struct menu { enum OutputFormat { PLAIN, SHELL, JSON }; +/* Footprint of input row spans for the output map, built in footprint.c. The + * struct is private to footprint.c. */ +struct footprint; +extern struct footprint * +fp_create(const struct Cell_head *ohd, const struct Cell_head *ihd, + const struct pj_info *oproj, const struct pj_info *iproj, + const struct pj_info *tproj, const double *y_center); +extern void fp_span(const struct footprint *g, int first_row, int end_row, + int first_col, int end_col, int *imin, int *imax); +extern int fp_num_blocks(const struct footprint *g); +extern int fp_block_start(const struct footprint *g, int b); +extern int fp_band_geometry(const struct footprint *g, int first_row, + size_t cap_bytes, int out_mult, int cell_size, + int in_cols, int *tile_blocks_out, + int *worst_block_rows); +extern void fp_free(struct footprint *g); + extern void bordwalk(const struct Cell_head *, struct Cell_head *, const struct pj_info *, const struct pj_info *, const struct pj_info *, int); @@ -67,7 +99,19 @@ extern void p_lanczos(struct cache *, void *, int, double, double, extern void p_lanczos_f(struct cache *, void *, int, double, double, struct Cell_head *); -#if 1 +/* interp_strip.c - strip versions of the resampling methods */ +extern void strip_bilinear(const struct strip *, void *, int, double, double, + struct Cell_head *, int *, int *); +extern void strip_cubic(const struct strip *, void *, int, double, double, + struct Cell_head *, int *, int *); +extern void strip_lanczos(const struct strip *, void *, int, double, double, + struct Cell_head *, int *, int *); +extern void strip_bilinear_f(const struct strip *, void *, int, double, double, + struct Cell_head *, int *, int *); +extern void strip_cubic_f(const struct strip *, void *, int, double, double, + struct Cell_head *, int *, int *); +extern void strip_lanczos_f(const struct strip *, void *, int, double, double, + struct Cell_head *, int *, int *); #define BKIDX(c, y, x) ((y) * (c)->stride + (x)) #define BKPTR(c, y, x) ((c)->grid[BKIDX((c), (y), (x))]) @@ -76,28 +120,4 @@ extern void p_lanczos_f(struct cache *, void *, int, double, double, : get_block((c), BKIDX((c), (y), (x)))) #define CVAL(c, y, x) ((*BLOCK((c), HI((y)), HI((x))))[LO((y))][LO((x))]) -#else - -static inline int BKIDX(const struct cache *c, int y, int x) -{ - return y * c->stride + x; -} - -static inline block *BKPTR(const struct cache *c, int y, int x) -{ - return c->grid[BKIDX(c, y, x)]; -} - -static inline block *BLOCK(struct cache *c, int y, int x) -{ - return BKPTR(c, y, x) ? BKPTR(c, y, x) : get_block(c, BKIDX(c, y, x)); -} - -static inline FCELL *CPTR(struct cache *c, int y, int x) -{ - return &(*BLOCK(c, HI(y), HI(x)))[LO(y)][LO(x)]; -} - -#endif - #endif diff --git a/raster/r.proj/r.proj.html b/raster/r.proj/r.proj.html index 80bf9f8efa3..b5bf020ad27 100644 --- a/raster/r.proj/r.proj.html +++ b/raster/r.proj/r.proj.html @@ -188,10 +188,19 @@

NOTES

r.proj is then run for the raster map the user wants to reproject. In this case a little preparation goes a long way.

+The nprocs parameter sets the number of compute threads +(default 1) and the memory parameter sets a limit in MB +(default 300). The module reads the input in bands sized to fit the +memory option, and threads split the output rows. A band is read +as column tiles when a full row does not fit. When nothing fits it falls +back to the serial cache with a message saying how much memory the map +needs. The sizing table itself is small and is not counted against +memory. To use the parallel path GRASS must be compiled with +OpenMP. +

When reprojecting whole-world maps the user should disable map-trimming with the -n flag. Trimming is not useful here -because the module has the whole map in memory anyway. Besides that, -world "edges" are hard (or impossible) to find in CRSs other +because world "edges" are hard (or impossible) to find in CRSs other than latitude-longitude so results may be odd with trimming.

EXAMPLES

diff --git a/raster/r.proj/r.proj.md b/raster/r.proj/r.proj.md index ece211a3902..79091383711 100644 --- a/raster/r.proj/r.proj.md +++ b/raster/r.proj/r.proj.md @@ -171,11 +171,19 @@ geodetic length of a pixel). *r.proj* is then run for the raster map the user wants to reproject. In this case a little preparation goes a long way. +The **nprocs** parameter sets the number of compute threads (default 1) +and the **memory** parameter sets a limit in MB (default 300). The module +reads the input in bands sized to fit the **memory** option, and threads +split the output rows. A band is read as column tiles when a full row does +not fit. When nothing fits it falls back to the serial cache with a message +saying how much memory the map needs. The sizing table itself is small and +is not counted against **memory**. To use the parallel path GRASS must be +compiled with OpenMP. + When reprojecting whole-world maps the user should disable map-trimming -with the **-n** flag. Trimming is not useful here because the module has -the whole map in memory anyway. Besides that, world "edges" are hard (or -impossible) to find in CRSs other than latitude-longitude so results may -be odd with trimming. +with the **-n** flag. Trimming is not useful here because world "edges" +are hard (or impossible) to find in CRSs other than latitude-longitude so +results may be odd with trimming. ## EXAMPLES diff --git a/raster/r.proj/tests/conftest.py b/raster/r.proj/tests/conftest.py index f7813f2a6d8..eee96685a31 100644 --- a/raster/r.proj/tests/conftest.py +++ b/raster/r.proj/tests/conftest.py @@ -41,3 +41,13 @@ def session_3857(gisdbase_with_source): gisdbase_with_source / "dst3857", env=os.environ.copy() ) as session: yield session + + +@pytest.fixture(scope="session") +def session_pole(gisdbase_with_source): + """Active session in an EPSG:3413 (north polar stereographic) project.""" + gs.create_project(gisdbase_with_source / "dst_pole", epsg="3413") + with gs.setup.init( + gisdbase_with_source / "dst_pole", env=os.environ.copy() + ) as session: + yield session diff --git a/raster/r.proj/tests/r_proj_parallel_test.py b/raster/r.proj/tests/r_proj_parallel_test.py new file mode 100644 index 00000000000..a7c38503d1c --- /dev/null +++ b/raster/r.proj/tests/r_proj_parallel_test.py @@ -0,0 +1,221 @@ +"""Parallel correctness tests for r.proj. Each test compares the +module's own nprocs=1 run against a multithreaded run, and the fallback +test forces the tile cache path.""" + +import pytest + +import grass.script as gs + +# Mirror of the names created in conftest.py. +SRC_PROJECT = "src4326" +INPUT_MID = "input_mid" +INPUT_POLAR = "input_polar" + + +def _env(session, **overrides): + """Session env copy with per-run overrides; never mutates the original.""" + env = dict(session.env) + for key, value in overrides.items(): + env[key] = str(value) + return env + + +# The suggested bounds depend only on the input, not the method, so cache them +# per input and avoid an r.proj -g call in every test. +_region_cache = {} + + +def _set_region_from_source(env, input_raster, method): + """Set the output region to r.proj's suggested bounds for the input.""" + if input_raster not in _region_cache: + text = gs.read_command( + "r.proj", + project=SRC_PROJECT, + mapset="PERMANENT", + input=input_raster, + method=method, + flags="g", + env=env, + ) + _region_cache[input_raster] = dict(token.split("=") for token in text.split()) + region = _region_cache[input_raster] + gs.run_command( + "g.region", + n=region["n"], + s=region["s"], + e=region["e"], + w=region["w"], + rows=region["rows"], + cols=region["cols"], + env=env, + ) + + +def _project(env, input_raster, output, method, **extra): + gs.run_command( + "r.proj", + project=SRC_PROJECT, + mapset="PERMANENT", + input=input_raster, + output=output, + method=method, + overwrite=True, + quiet=True, + env=env, + **extra, + ) + + +def _stats(env, raster): + return gs.parse_command("r.univar", map=raster, flags="g", env=env) + + +def _assert_bitwise_identical(env, a, b, diff): + """Check a and b are bitwise identical and null in the same cells.""" + gs.run_command( + "r.mapcalc", expression=f"{diff} = abs({a} - {b})", overwrite=True, env=env + ) + gs.run_command( + "r.mapcalc", + expression=f"{diff}_null = if(isnull({a}) != isnull({b}), 1, 0)", + overwrite=True, + env=env, + ) + sa = _stats(env, a) + sb = _stats(env, b) + sd = _stats(env, diff) + assert int(sa["n"]) > 0, "output is empty; the comparison would be vacuous" + assert int(sa["n"]) == int(sb["n"]) + assert int(sa["null_cells"]) == int(sb["null_cells"]) + assert float(sd["max"]) == 0.0 + # abs(a - b) is null wherever either map is null, so a swapped null hides + # from sd. The null map catches a null that moved. + assert float(_stats(env, f"{diff}_null")["max"]) == 0.0 + + +def test_bilinear_parallel_matches_serial(session_3857): + """Bilinear method needs to match serial bitwise. It first checks that + bilinear and nearest outputs differ, so a silent fallback to nearest + cannot happen.""" + session = session_3857 + base = _env(session) + _set_region_from_source(base, INPUT_MID, "bilinear") + + _project(base, INPUT_MID, "bilin_serial", "bilinear", nprocs=1) + _project(base, INPUT_MID, "nearest_ref", "nearest", nprocs=1) + gs.run_command( + "r.mapcalc", + expression="dispatch_live = abs(bilin_serial - nearest_ref)", + overwrite=True, + env=base, + ) + assert float(_stats(base, "dispatch_live")["max"]) > 0, ( + "bilinear output equals nearest; dispatch may have fallen back" + ) + + _project(base, INPUT_MID, "bilin_parallel", "bilinear", nprocs=4) + _assert_bitwise_identical(base, "bilin_serial", "bilin_parallel", "bilin_diff") + + +def test_nearest_low_memory_matches_serial(session_3857): + """The nearest method at a small memory cap (memory=5, nprocs=4) matches the + default memory serial run bitwise. The 50x50 fixture fits one band even at + this cap, so this checks the low-memory path rather than multi-band.""" + session = session_3857 + base = _env(session) + _set_region_from_source(base, INPUT_MID, "nearest") + + _project(base, INPUT_MID, "mem_serial", "nearest", nprocs=1) + _project(base, INPUT_MID, "mem_low", "nearest", nprocs=4, memory=5) + _assert_bitwise_identical(base, "mem_serial", "mem_low", "mem_diff") + + +def test_pole_nearest_parallel_matches_serial(session_pole): + """Nearest at the north pole matches serial bitwise.""" + session = session_pole + base = _env(session) + # Fixed 1200 km box centered on the pole (EPSG:3413 meters), 50x50. + gs.run_command( + "g.region", + n=600000, + s=-600000, + e=600000, + w=-600000, + rows=50, + cols=50, + env=base, + ) + + _project(base, INPUT_POLAR, "pole_serial", "nearest", nprocs=1) + _project(base, INPUT_POLAR, "pole_parallel", "nearest", nprocs=4) + _assert_bitwise_identical(base, "pole_serial", "pole_parallel", "pole_diff") + + +@pytest.mark.parametrize("method", ["nearest", "bilinear", "lanczos"]) +def test_forced_fallback_matches_banded(session_3857, method): + """Forcing the tile cache with R_PROJ_FORCE_TILECACHE=1 gives the same output + as the banded path. One method per kernel family keeps the suite fast. This + pairs the legacy p_ cache kernels against the strip kernels.""" + session = session_3857 + base = _env(session) + _set_region_from_source(base, INPUT_MID, method) + + _project( + _env(session, R_PROJ_FORCE_TILECACHE=1), + INPUT_MID, + "fallback_tilecache", + method, + nprocs=1, + ) + _project(base, INPUT_MID, "banded", method, nprocs=4) + _assert_bitwise_identical(base, "fallback_tilecache", "banded", "fallback_diff") + + +def _project_capture(env, input_raster, output, method, **extra): + """Run r.proj and return the messages it writes to stderr.""" + env = dict(env, GRASS_VERBOSE="3", LC_ALL="C") + proc = gs.start_command( + "r.proj", + project=SRC_PROJECT, + mapset="PERMANENT", + input=input_raster, + output=output, + method=method, + overwrite=True, + env=env, + stderr=gs.PIPE, + **extra, + ) + stderr = proc.communicate()[1] + if isinstance(stderr, bytes): + stderr = stderr.decode() + assert proc.returncode == 0, stderr + return stderr + + +def test_forced_reload_matches_serial(session_3857): + """FG_SHRINK_SPAN shrinks each band so the reload runs, and the reloaded run + still matches the nprocs=1 run bitwise.""" + session = session_3857 + base = _env(session) + _set_region_from_source(base, INPUT_MID, "nearest") + + _project(base, INPUT_MID, "reload_serial", "nearest", nprocs=1) + + shrunk = _project_capture( + _env(session, FG_SHRINK_SPAN=2), + INPUT_MID, + "reload_shrunk", + "nearest", + nprocs=4, + ) + # The shrunk strips leave out rows the projection needs, so r.proj reloads + # and logs it here. Without this marker the run never hit the reload path. + assert "Reloading band strip" in shrunk + + # The full-size strips already cover every needed row, so this run never + # reloads. That is what ties the marker above to the shrink. + unshrunk = _project_capture(base, INPUT_MID, "reload_unshrunk", "nearest", nprocs=4) + assert "Reloading band strip" not in unshrunk + + _assert_bitwise_identical(base, "reload_serial", "reload_shrunk", "reload_diff")