diff --git a/raster/CMakeLists.txt b/raster/CMakeLists.txt index cb8a741043d..ad7bdeb8611 100644 --- a/raster/CMakeLists.txt +++ b/raster/CMakeLists.txt @@ -244,6 +244,7 @@ build_program_in_subdir( build_program_in_subdir( r.geomorphon DEPENDS ${LIBM} grass_gis grass_gmath grass_parson grass_raster + OPTIONAL_DEPENDS OpenMP::OpenMP_C ) build_program_in_subdir(r.grow.distance DEPENDS grass_gis grass_raster ${LIBM}) diff --git a/raster/r.geomorphon/Makefile b/raster/r.geomorphon/Makefile index 8845a0297fb..1ed4890bfbb 100644 --- a/raster/r.geomorphon/Makefile +++ b/raster/r.geomorphon/Makefile @@ -3,7 +3,10 @@ MODULE_TOPDIR = ../.. PGM = r.geomorphon LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB) +EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB) DEPENDENCIES = $(RASTERDEP) $(GISDEP) +EXTRA_CFLAGS = $(OPENMP_CFLAGS) +EXTRA_INC = $(OPENMP_INCPATH) include $(MODULE_TOPDIR)/include/Make/Module.make diff --git a/raster/r.geomorphon/benchmark/benchmark_r_geomorphon.py b/raster/r.geomorphon/benchmark/benchmark_r_geomorphon.py new file mode 100644 index 00000000000..bbc5b0a57fa --- /dev/null +++ b/raster/r.geomorphon/benchmark/benchmark_r_geomorphon.py @@ -0,0 +1,96 @@ +"""Benchmarking of r.geomorphon +raster (2D) +""" + +from grass.exceptions import CalledModuleError, GrassError +from grass.pygrass.modules import Module +import grass.benchmark as bm + +# Baselines held fixed while one dimension is swept. +BASE_MAPSIZE = 50e6 # cells +BASE_SEARCH = 25 # outer search radius in cells +BASE_MEMORY = 300 # MB +MAPSIZES = [10e6, 50e6, 100e6] +SEARCHES = [10, 25, 50] +METRICS = ["time", "speedup", "efficiency"] + + +def main(): + # Sweep raster size at the baseline search radius. + results = [] + for mapsize in MAPSIZES: + benchmark( + size=int(mapsize**0.5), + search=BASE_SEARCH, + memory=BASE_MEMORY, + label=f"r.geomorphon_{int(mapsize / 1e6)}M", + results=results, + ) + plot(results, "rastersize") + + # Sweep search radius at the baseline raster size. + results = [] + for search in SEARCHES: + benchmark( + size=int(BASE_MAPSIZE**0.5), + search=search, + memory=BASE_MEMORY, + label=f"r.geomorphon_search_{search}", + results=results, + ) + plot(results, "search") + + +def benchmark(size, search, memory, label, results): + reference = "benchmark_r_geomorphon_reference" + output = "benchmark_r_geomorphon" + generate_map(rows=size, cols=size, fname=reference) + module = Module( + "r.geomorphon", + elevation=reference, + forms=output, + search=search, + memory=memory, + run_=False, + overwrite=True, + ) + results.append( + bm.benchmark_nprocs( + module, + label=label, + max_nprocs=8, + repeat=3, + ) + ) + Module( + "g.remove", + quiet=True, + flags="f", + type="raster", + pattern="benchmark_r_geomorphon*", + ) + + +def plot(results, sweep): + for metric in METRICS: + bm.nprocs_plot( + results, + filename=f"r_geomorphon_{sweep}_{metric}.svg", + title=f"r.geomorphon {sweep} {metric}", + metric=metric, + ) + + +def generate_map(rows, cols, fname): + Module("g.region", flags="p", n=rows, e=cols, res=1, w=0, s=0) + # Generate using r.random.surface if r.surf.fractal fails + try: + print("Generating reference map using r.surf.fractal...") + Module("r.surf.fractal", output=fname, overwrite=True) + except (CalledModuleError, GrassError): + print("r.surf.fractal fails, using r.random.surface instead...") + Module("r.random.surface", output=fname, overwrite=True) + + +if __name__ == "__main__": + main() diff --git a/raster/r.geomorphon/local_proto.h b/raster/r.geomorphon/local_proto.h index 5cfe0e31c6b..aa4267240f6 100644 --- a/raster/r.geomorphon/local_proto.h +++ b/raster/r.geomorphon/local_proto.h @@ -32,7 +32,6 @@ typedef struct { char elevname[150]; RASTER_MAP_TYPE raster_type; - FCELL **elev; int fd; /* file descriptor */ } MAPS; @@ -64,11 +63,24 @@ typedef enum { CNT /* counter */ } FORMS; +/* Invariant compute_forms inputs shared by every cell in a run. */ +struct geomorphon_config { + double search_dist, skip_dist, flat_dist, max_resolution; + int extended, oneoff; +}; + +/* Per-cell compute_forms outputs. */ +struct geomorphon_result { + PATTERN *pattern; + int pattern_size; + FORMS cur_form, orig_form; + double eff_search, eff_skip, eff_flat; +}; + /* main */ GLOBAL MAPS elevation; GLOBAL int ncols, row_radius_size, row_buffer_size; GLOBAL int skip_cells; -GLOBAL double search_distance, flat_distance; GLOBAL double flat_threshold, flat_threshold_height; GLOBAL struct Cell_head window; @@ -77,15 +89,16 @@ GLOBAL enum { ANGLEV1, ANGLEV2, ANGLEV2_DISTANCE } compmode; /* memory */ int open_map(MAPS *rast); -int shift_buffers(int row); -int free_map(FCELL **map, int n); +int load_strip(int fd, RASTER_MAP_TYPE rtype, void *tmp_buf, FCELL **rows, + int abs_first, int count); int write_form_cat_colors(char *raster); int write_contrast_colors(char *); const char *form_short_name(const FORMS); const char *form_long_name(const FORMS); /* pattern */ -int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, const int); +int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, const int, + double search_distance, double flat_distance, FCELL **rows); extern const char *direction_name[]; /* geom */ diff --git a/raster/r.geomorphon/main.c b/raster/r.geomorphon/main.c index b21a62a07ac..0ae1977c26e 100644 --- a/raster/r.geomorphon/main.c +++ b/raster/r.geomorphon/main.c @@ -21,13 +21,15 @@ * *****************************************************************************/ +#if defined(_OPENMP) +#include +#endif + #include #include #define MAIN #include "local_proto.h" -#define WINDOW_THRESHOLD 100000000 - typedef enum { o_forms, o_ternary, @@ -44,6 +46,12 @@ typedef enum { o_size } outputs; +/* Compute one cell's pattern and landform, shared by the raster and one-off + * paths. */ +static void compute_forms(FCELL **rows, int cur_row, int row, int col, + const struct geomorphon_config *cfg, + PATTERN *patterns, struct geomorphon_result *res); + int main(int argc, char **argv) { struct { /* struct is used both for interface and output */ @@ -87,7 +95,8 @@ int main(int argc, char **argv) struct GModule *module; struct Option *opt_input, *opt_output[o_size], *par_search_radius, *par_skip_radius, *par_flat_threshold, *par_flat_distance, - *par_comparison, *par_coords, *par_profiledata, *par_profileformat; + *par_comparison, *par_coords, *par_profiledata, *par_profileformat, + *par_memory, *par_nprocs; struct Flag *flag_units, *flag_extended; struct History history; @@ -95,10 +104,12 @@ int main(int argc, char **argv) int i; int meters = 0, extended = 0; /* flags */ int oneoff; - int row, cur_row, col; int nrows; - int pattern_size; + int brows; /* rows per output band buffer */ + int memory; /* memory cap in MB from the memory= option */ + int nprocs; /* number of OpenMP threads */ int search_cells; + double search_distance, flat_distance; double skip_distance; double max_resolution; double oneoff_easting, oneoff_northing; @@ -198,6 +209,9 @@ int main(int argc, char **argv) G_option_requires(par_profileformat, par_coords, NULL); G_option_requires(par_coords, par_profileformat, NULL); + par_nprocs = G_define_standard_option(G_OPT_M_NPROCS); + par_memory = G_define_standard_option(G_OPT_MEMORYMB); + if (G_parser(argc, argv)) exit(EXIT_FAILURE); } @@ -290,6 +304,42 @@ int main(int argc, char **argv) "At least %d rows are needed. Set larger computational " "region with g.region or use a smaller search value."), par_search_radius->answer, nrows, row_buffer_size + 1); + /* Band height from the memory cap. */ + memory = atoi(par_memory->answer); + nprocs = G_set_omp_num_threads(par_nprocs); + nprocs = Rast_disable_omp_on_mask(nprocs); + { + size_t fixed = (size_t)2 * row_radius_size * ncols * sizeof(FCELL); + size_t cap = (size_t)memory * (1 << 20); + size_t per_row = (size_t)ncols * sizeof(FCELL) * (1 + num_outputs); + size_t budget; + + size_t need_bytes = fixed + (size_t)ncols * sizeof(FCELL) + + (size_t)nprocs * per_row; + if (cap < need_bytes) { + long need = (long)((need_bytes + (1 << 20) - 1) >> 20); + + G_warning( + _("memory=%d MB is below the %ld MB minimum for " + "search=%s. This run will finish with the minimum MB " + "needed instead."), + memory, need, par_search_radius->answer); + budget = (cap > fixed) ? cap - fixed : 0; + } + else + budget = cap - fixed; + size_t brows_sz = budget / per_row; + + if (brows_sz > (size_t)nrows) + brows = nrows; + else + brows = (int)brows_sz; + if (brows < nprocs) + brows = nprocs; + } + /* One-off needs only its own row, so a single band suffices. */ + if (oneoff) + brows = 1; search_distance = (meters) ? search_radius : ns_resolution * search_cells; @@ -333,18 +383,6 @@ int main(int argc, char **argv) * to the cell boundary). And also perhaps another "-r" flag to * restore the region afterwards. */ - if (oneoff) { - unsigned long window_square = nrows * ncols; - unsigned long search_square = 4 * search_cells * search_cells; - - if (window_square > WINDOW_THRESHOLD && - window_square / search_square > 10) - G_warning( - _("There may be a notable processing delay because the " - "computational region is %lu times larger than " - "necessary"), - window_square / search_square); - } } generate_ternary_codes(); @@ -353,285 +391,372 @@ int main(int argc, char **argv) strcpy(elevation.elevname, opt_input->answer); open_map(&elevation); - if (1) { - PATTERN *pattern; - PATTERN patterns[4]; - void *pointer_buf; + { double search_dist = search_distance; double skip_dist = skip_distance; double flat_dist = flat_distance; + struct geomorphon_config cfg = {search_dist, skip_dist, flat_dist, + max_resolution, extended, oneoff}; double area_of_octagon = 4 * (search_distance * search_distance) * sin(DEGREE2RAD(45.)); unsigned char oneoff_done = 0; + /* Shared input strip, reused for every band. */ + int strip_cap = brows + 2 * row_radius_size + 1; + FCELL *strip_block = + (FCELL *)G_malloc((size_t)strip_cap * ncols * sizeof(FCELL)); + FCELL **strip_ptr = + (FCELL **)G_malloc((size_t)strip_cap * sizeof(FCELL *)); - /* prepare outputs */ - for (i = o_forms; i < o_size; ++i) - if (opt_output[i]->answer) { - rasters[i].fd = Rast_open_new(opt_output[i]->answer, - rasters[i].out_data_type); - rasters[i].buffer = Rast_allocate_buf(rasters[i].out_data_type); - } - - /* main loop */ - for (row = 0; row < nrows; ++row) { - G_percent(row, nrows, 2); - cur_row = (row < row_radius_size) - ? row - : ((row >= nrows - row_radius_size - 1) - ? row_buffer_size - (nrows - row - 1) - : row_radius_size); - - if (row > (row_radius_size) && row < nrows - (row_radius_size + 1)) - shift_buffers(row); - - /* If skipping the current row, only after the buffer shift would be - * fine. */ - if (oneoff && row != oneoff_row) - continue; - - for (col = 0; col < ncols; ++col) { - /* If skipping the current column, early would be fine. */ - if (oneoff && col != oneoff_col) - continue; - - /* on borders forms ussualy are innatural. */ - if (row < (skip_cells + 1) || row > nrows - (skip_cells + 2) || - col < (skip_cells + 1) || col > ncols - (skip_cells + 2) || - Rast_is_f_null_value(&elevation.elev[cur_row][col])) { - /* set outputs to NULL and do nothing if source value is - * null or border */ - for (i = o_forms; i < o_size; ++i) - if (opt_output[i]->answer) { - pointer_buf = rasters[i].buffer; - switch (rasters[i].out_data_type) { - case CELL_TYPE: - Rast_set_c_null_value( - &((CELL *)pointer_buf)[col], 1); - break; - case FCELL_TYPE: - Rast_set_f_null_value( - &((FCELL *)pointer_buf)[col], 1); - break; - case DCELL_TYPE: - Rast_set_d_null_value( - &((DCELL *)pointer_buf)[col], 1); - break; - default: - G_fatal_error(_("Unknown output data type")); - } - } - continue; - } /* end null value */ - { - FORMS cur_form; - FORMS orig_form; - - search_distance = search_dist; - skip_distance = skip_dist; - flat_distance = flat_dist; - pattern_size = - calc_pattern(&patterns[0], row, cur_row, col, oneoff); - pattern = &patterns[0]; - cur_form = orig_form = determine_form( - pattern->num_negatives, pattern->num_positives); - - /* correction of forms */ - if (extended && search_distance > 10 * max_resolution) { - /* 1) remove extensive innatural forms: ridges, peaks, - * shoulders and footslopes */ - if ((cur_form == SH || cur_form == FS || - cur_form == PK || cur_form == RI)) { - FORMS small_form; - - search_distance = - (search_dist / 2. < 4 * max_resolution) - ? 4 * max_resolution - : search_dist / 2.; - skip_distance = 0; - flat_distance = 0; - pattern_size = calc_pattern(&patterns[1], row, - cur_row, col, 0); - pattern = &patterns[1]; - small_form = determine_form(pattern->num_negatives, - pattern->num_positives); - if (cur_form == SH || cur_form == FS) - cur_form = (small_form == FL) ? FL : cur_form; - if (cur_form == PK || cur_form == RI) - cur_form = small_form; - } - /* 3) Depressions */ - - } /* end of correction */ - - /* one-off mode */ - if (oneoff) { - char buf[BUFSIZ]; - float azimuth, elongation, width; - - radial2cartesian(pattern); - shape(pattern, pattern_size, &azimuth, &elongation, - &width); - prof_map_info(); - prof_sso("computation_parameters"); - prof_dbl("easting", oneoff_easting); - prof_dbl("northing", oneoff_northing); - prof_mtr("search_m", search_distance); - prof_int("search_cells", search_cells); - prof_mtr("skip_m", skip_distance); - prof_int("skip_cells", skip_cells); - prof_dbl("flat_thresh_deg", RAD2DEGREE(flat_threshold)); - prof_mtr("flat_distance_m", flat_distance); - prof_mtr("flat_height_m", flat_threshold_height); - prof_bln("extended_correction", extended); - prof_eso(); /* computation_parameters */ - prof_sso("intermediate_data"); - if (extended) { - prof_int("initial_landform_cat", orig_form); - prof_str("initial_landform_code", - form_short_name(orig_form)); - prof_str("initial_landform_name", - form_long_name(orig_form)); - } - prof_int("ternary_498", - determine_ternary(pattern->pattern)); - prof_int("ternary_6561", - preliminary_ternary(pattern->pattern)); - prof_int("pattern_size", pattern_size); - prof_dbl("origin_easting", - Rast_col_to_easting(col + 0.5, &window)); - prof_dbl("origin_northing", - Rast_row_to_northing(row + 0.5, &window)); - prof_pattern(elevation.elev[cur_row][col], pattern); - prof_eso(); /* intermediate_data */ - prof_sso("final_results"); - prof_int("landform_cat", cur_form); - prof_str("landform_code", form_short_name(cur_form)); - prof_str("landform_name", form_long_name(cur_form)); - prof_int("landform_deviation", - form_deviation(pattern->num_negatives, - pattern->num_positives)); - prof_dbl("azimuth", azimuth); - prof_dbl("elongation", elongation); - prof_mtr("width_m", width); - prof_mtr("intensity_m", - intensity(pattern->elevation, pattern_size)); - prof_mtr("exposition_m", - exposition(pattern->elevation)); - prof_mtr("range_m", range(pattern->elevation)); - prof_dbl("variance", - variance(pattern->elevation, pattern_size)); - prof_dbl("extends", extends(pattern) / area_of_octagon); - prof_mtr("octagon_perimeter_m", - octa_perimeter(pattern)); - prof_mtr("octagon_area_m2", extends(pattern)); - prof_mtr("mesh_perimeter_m", mesh_perimeter(pattern)); - prof_mtr("mesh_area_m2", mesh_area(pattern)); - prof_eso(); /* final_results */ - /* - * When adding new data items, increment the minor - * version. When deleting, moving, renaming or - * otherwise changing existing data, increment the - * major version and reset the minor version. - */ - prof_int("format_version_major", 0); - prof_int("format_version_minor", 9); - prof_utc("timestamp", time(NULL)); - snprintf(buf, sizeof(buf), "r.geomorphon GRASS %s [%s]", - GRASS_VERSION_STRING, GRASS_HEADERS_VERSION); - prof_str("generator", buf); - - oneoff_done = - prof_write(profile_file, par_profileformat->answer); - if (oneoff_done) - G_verbose_message( - _("Profile data has been written")); - else - G_important_message( - _("Failed writing profile data")); - /* Break out of both loops. */ - row = nrows - 1; - break; - } /* end of one-off mode */ - - pattern = &patterns[0]; - if (opt_output[o_forms]->answer) - ((CELL *)rasters[o_forms].buffer)[col] = cur_form; + if (oneoff) { + /* Serial single-cell profile, outside any parallel region. */ + void *tmp_buf = Rast_allocate_buf(elevation.raster_type); + int win_top = MIN(MAX(oneoff_row - row_radius_size, 0), + nrows - row_buffer_size - 1); + int strip_lo = win_top; + int strip_hi = win_top + row_buffer_size; + int strip_rows, cur_row = oneoff_row - win_top, k; + FCELL **rows; + + if (strip_hi > nrows - 1) + strip_hi = nrows - 1; + strip_rows = strip_hi - strip_lo + 1; + for (k = 0; k < strip_rows; ++k) + strip_ptr[k] = strip_block + (size_t)k * ncols; + load_strip(elevation.fd, elevation.raster_type, tmp_buf, strip_ptr, + strip_lo, strip_rows); + rows = &strip_ptr[win_top - strip_lo]; + + if (!(oneoff_row < skip_cells + 1 || + oneoff_row > nrows - (skip_cells + 2) || + oneoff_col < skip_cells + 1 || + oneoff_col > ncols - (skip_cells + 2) || + Rast_is_f_null_value(&rows[cur_row][oneoff_col]))) { + PATTERN patterns[4]; + struct geomorphon_result res; + char buf[BUFSIZ]; + float azimuth, elongation, width; + + compute_forms(rows, cur_row, oneoff_row, oneoff_col, &cfg, + patterns, &res); + PATTERN *pattern = res.pattern; + int pattern_size = res.pattern_size; + FORMS cur_form = res.cur_form, orig_form = res.orig_form; + double eff_search = res.eff_search, eff_skip = res.eff_skip, + eff_flat = res.eff_flat; + radial2cartesian(pattern); + shape(pattern, pattern_size, &azimuth, &elongation, &width); + prof_map_info(); + prof_sso("computation_parameters"); + prof_dbl("easting", oneoff_easting); + prof_dbl("northing", oneoff_northing); + prof_mtr("search_m", eff_search); + prof_int("search_cells", search_cells); + prof_mtr("skip_m", eff_skip); + prof_int("skip_cells", skip_cells); + prof_dbl("flat_thresh_deg", RAD2DEGREE(flat_threshold)); + prof_mtr("flat_distance_m", eff_flat); + prof_mtr("flat_height_m", flat_threshold_height); + prof_bln("extended_correction", extended); + prof_eso(); /* computation_parameters */ + prof_sso("intermediate_data"); + if (extended) { + prof_int("initial_landform_cat", orig_form); + prof_str("initial_landform_code", + form_short_name(orig_form)); + prof_str("initial_landform_name", + form_long_name(orig_form)); } - - if (opt_output[o_ternary]->answer) - ((CELL *)rasters[o_ternary].buffer)[col] = - determine_ternary(pattern->pattern); - if (opt_output[o_positive]->answer) - ((CELL *)rasters[o_positive].buffer)[col] = - rotate(pattern->positives); - if (opt_output[o_negative]->answer) - ((CELL *)rasters[o_negative].buffer)[col] = - rotate(pattern->negatives); - if (opt_output[o_intensity]->answer) - ((FCELL *)rasters[o_intensity].buffer)[col] = - intensity(pattern->elevation, pattern_size); - if (opt_output[o_exposition]->answer) - ((FCELL *)rasters[o_exposition].buffer)[col] = - exposition(pattern->elevation); - if (opt_output[o_range]->answer) - ((FCELL *)rasters[o_range].buffer)[col] = - range(pattern->elevation); - if (opt_output[o_variance]->answer) - ((FCELL *)rasters[o_variance].buffer)[col] = - variance(pattern->elevation, pattern_size); - - /* used only for next four shape functions + prof_int("ternary_498", determine_ternary(pattern->pattern)); + prof_int("ternary_6561", preliminary_ternary(pattern->pattern)); + prof_int("pattern_size", pattern_size); + prof_dbl("origin_easting", + Rast_col_to_easting(oneoff_col + 0.5, &window)); + prof_dbl("origin_northing", + Rast_row_to_northing(oneoff_row + 0.5, &window)); + prof_pattern(rows[cur_row][oneoff_col], pattern); + prof_eso(); /* intermediate_data */ + prof_sso("final_results"); + prof_int("landform_cat", cur_form); + prof_str("landform_code", form_short_name(cur_form)); + prof_str("landform_name", form_long_name(cur_form)); + prof_int("landform_deviation", + form_deviation(pattern->num_negatives, + pattern->num_positives)); + prof_dbl("azimuth", azimuth); + prof_dbl("elongation", elongation); + prof_mtr("width_m", width); + prof_mtr("intensity_m", + intensity(pattern->elevation, pattern_size)); + prof_mtr("exposition_m", exposition(pattern->elevation)); + prof_mtr("range_m", range(pattern->elevation)); + prof_dbl("variance", + variance(pattern->elevation, pattern_size)); + prof_dbl("extends", extends(pattern) / area_of_octagon); + prof_mtr("octagon_perimeter_m", octa_perimeter(pattern)); + prof_mtr("octagon_area_m2", extends(pattern)); + prof_mtr("mesh_perimeter_m", mesh_perimeter(pattern)); + prof_mtr("mesh_area_m2", mesh_area(pattern)); + prof_eso(); /* final_results */ + /* + * When adding new data items, increment the minor version. + * When deleting, moving, renaming or otherwise changing + * existing data, increment the major version and reset the + * minor version. */ - if (opt_output[o_elongation]->answer || - opt_output[o_azimuth]->answer || - opt_output[o_extend]->answer || - opt_output[o_width]->answer) { - float azimuth, elongation, width; - - radial2cartesian(pattern); - shape(pattern, pattern_size, &azimuth, &elongation, &width); - if (opt_output[o_azimuth]->answer) - ((FCELL *)rasters[o_azimuth].buffer)[col] = azimuth; - if (opt_output[o_elongation]->answer) - ((FCELL *)rasters[o_elongation].buffer)[col] = - elongation; - if (opt_output[o_width]->answer) - ((FCELL *)rasters[o_width].buffer)[col] = width; + prof_int("format_version_major", 0); + prof_int("format_version_minor", 9); + prof_utc("timestamp", time(NULL)); + snprintf(buf, sizeof(buf), "r.geomorphon GRASS %s [%s]", + GRASS_VERSION_STRING, GRASS_HEADERS_VERSION); + prof_str("generator", buf); + oneoff_done = + prof_write(profile_file, par_profileformat->answer); + if (oneoff_done) + G_verbose_message(_("Profile data has been written")); + else + G_important_message(_("Failed writing profile data")); + } + G_free(tmp_buf); + } + else { + /* Parallel raster path, one shared strip loaded by all threads. */ + int *fd_thread = G_malloc(sizeof(int) * nprocs); + void **tmp_thread = G_malloc(sizeof(void *) * nprocs); + int computed = 0; + int band_start, t; + + for (t = 0; t < nprocs; ++t) { + fd_thread[t] = (t == 0) ? elevation.fd + : Rast_open_old(elevation.elevname, ""); + tmp_thread[t] = Rast_allocate_buf(elevation.raster_type); + } + for (i = o_forms; i < o_size; ++i) + if (opt_output[i]->answer) { + rasters[i].fd = Rast_open_new(opt_output[i]->answer, + rasters[i].out_data_type); + rasters[i].buffer = + G_malloc((size_t)brows * ncols * + Rast_cell_size(rasters[i].out_data_type)); } - if (opt_output[o_extend]->answer) - ((FCELL *)rasters[o_extend].buffer)[col] = - extends(pattern) / area_of_octagon; - } /* end for col */ + for (band_start = 0; band_start < nrows; band_start += brows) { + int band_end = band_start + brows; + int band_count, strip_lo, strip_hi, strip_rows, k, r; + + if (band_end > nrows) + band_end = nrows; + band_count = band_end - band_start; + + strip_lo = MIN(MAX(band_start - row_radius_size, 0), + nrows - row_buffer_size - 1); + strip_hi = MIN(MAX(band_end - 1 - row_radius_size, 0), + nrows - row_buffer_size - 1) + + row_buffer_size; + if (strip_hi > nrows - 1) + strip_hi = nrows - 1; + strip_rows = strip_hi - strip_lo + 1; + for (k = 0; k < strip_rows; ++k) + strip_ptr[k] = strip_block + (size_t)k * ncols; + +#pragma omp parallel if (nprocs > 1) + { + int tid = 0, nth = 1; + int ls, le, rs, re, row; +#if defined(_OPENMP) + tid = omp_get_thread_num(); + nth = omp_get_num_threads(); +#endif + /* Each thread loads a disjoint chunk of the strip. */ + ls = strip_rows * tid / nth; + le = strip_rows * (tid + 1) / nth; + if (le > ls) + load_strip(fd_thread[tid], elevation.raster_type, + tmp_thread[tid], &strip_ptr[ls], + strip_lo + ls, le - ls); + + /* Barrier so all strip writes are visible before any thread + * reads. */ +#pragma omp barrier + + /* Split this band's output rows across threads. */ + rs = band_start + band_count * tid / nth; + re = band_start + band_count * (tid + 1) / nth; + for (row = rs; row < re; ++row) { + int win_top = MIN(MAX(row - row_radius_size, 0), + nrows - row_buffer_size - 1); + int cur_row = row - win_top; + FCELL **rows = &strip_ptr[win_top - strip_lo]; + size_t boff = (size_t)(row - band_start) * ncols; + int col; + + if (tid == 0) { + int done; +#pragma omp atomic read + done = computed; + G_percent(done, nrows, 2); + } - /* write existing outputs */ + for (col = 0; col < ncols; ++col) { + PATTERN patterns[4]; + PATTERN *pattern; + void *pointer_buf; + int pattern_size; + + /* on borders forms ussualy are innatural. */ + if (row < (skip_cells + 1) || + row > nrows - (skip_cells + 2) || + col < (skip_cells + 1) || + col > ncols - (skip_cells + 2) || + Rast_is_f_null_value(&rows[cur_row][col])) { + /* set outputs to NULL on null source or border + */ + int m; + + for (m = o_forms; m < o_size; ++m) + if (opt_output[m]->answer) { + pointer_buf = rasters[m].buffer; + switch (rasters[m].out_data_type) { + case CELL_TYPE: + Rast_set_c_null_value( + &((CELL *) + pointer_buf)[boff + col], + 1); + break; + case FCELL_TYPE: + Rast_set_f_null_value( + &((FCELL *) + pointer_buf)[boff + col], + 1); + break; + case DCELL_TYPE: + Rast_set_d_null_value( + &((DCELL *) + pointer_buf)[boff + col], + 1); + break; + default: + G_fatal_error( + _("Unknown output data type")); + } + } + continue; + } /* end null value */ + { + struct geomorphon_result res; + + compute_forms(rows, cur_row, row, col, &cfg, + patterns, &res); + pattern_size = res.pattern_size; + /* The extended correction only changes the + * forms output. Every other output is computed + * from the original full search pattern, + * exactly like the serial code does. */ + pattern = &patterns[0]; + if (opt_output[o_forms]->answer) + ((CELL *)rasters[o_forms] + .buffer)[boff + col] = res.cur_form; + if (opt_output[o_ternary]->answer) + ((CELL *)rasters[o_ternary] + .buffer)[boff + col] = + determine_ternary(pattern->pattern); + if (opt_output[o_positive]->answer) + ((CELL *)rasters[o_positive] + .buffer)[boff + col] = + rotate(pattern->positives); + if (opt_output[o_negative]->answer) + ((CELL *)rasters[o_negative] + .buffer)[boff + col] = + rotate(pattern->negatives); + if (opt_output[o_intensity]->answer) + ((FCELL *)rasters[o_intensity] + .buffer)[boff + col] = + intensity(pattern->elevation, + pattern_size); + if (opt_output[o_exposition]->answer) + ((FCELL *)rasters[o_exposition] + .buffer)[boff + col] = + exposition(pattern->elevation); + if (opt_output[o_range]->answer) + ((FCELL *)rasters[o_range] + .buffer)[boff + col] = + range(pattern->elevation); + if (opt_output[o_variance]->answer) + ((FCELL *)rasters[o_variance] + .buffer)[boff + col] = + variance(pattern->elevation, + pattern_size); + if (opt_output[o_elongation]->answer || + opt_output[o_azimuth]->answer || + opt_output[o_extend]->answer || + opt_output[o_width]->answer) { + float azimuth, elongation, width; + + radial2cartesian(pattern); + shape(pattern, pattern_size, &azimuth, + &elongation, &width); + if (opt_output[o_azimuth]->answer) + ((FCELL *)rasters[o_azimuth] + .buffer)[boff + col] = azimuth; + if (opt_output[o_elongation]->answer) + ((FCELL *)rasters[o_elongation] + .buffer)[boff + col] = elongation; + if (opt_output[o_width]->answer) + ((FCELL *)rasters[o_width] + .buffer)[boff + col] = width; + } + if (opt_output[o_extend]->answer) + ((FCELL *)rasters[o_extend] + .buffer)[boff + col] = + extends(pattern) / area_of_octagon; + } + } /* end for col */ +#pragma omp atomic update + computed++; + } /* end for row */ + } /* end parallel */ + + /* serial in-order drain of this band's rows, per map */ + for (i = o_forms; i < o_size; ++i) + if (opt_output[i]->answer) { + size_t rowsz = (size_t)ncols * + Rast_cell_size(rasters[i].out_data_type); + + for (r = 0; r < band_count; ++r) + Rast_put_row(rasters[i].fd, + (char *)rasters[i].buffer + + (size_t)r * rowsz, + rasters[i].out_data_type); + } + } /* end band loop */ + G_percent(nrows, nrows, 2); + + for (t = 0; t < nprocs; ++t) { + Rast_close(fd_thread[t]); + G_free(tmp_thread[t]); + } + G_free(fd_thread); + G_free(tmp_thread); for (i = o_forms; i < o_size; ++i) - if (opt_output[i]->answer) - Rast_put_row(rasters[i].fd, rasters[i].buffer, - rasters[i].out_data_type); + if (opt_output[i]->answer) { + G_free(rasters[i].buffer); + Rast_close(rasters[i].fd); + Rast_short_history(opt_output[i]->answer, "raster", + &history); + Rast_command_history(&history); + Rast_write_history(opt_output[i]->answer, &history); + } + if (opt_output[o_forms]->answer) + write_form_cat_colors(opt_output[o_forms]->answer); + if (opt_output[o_intensity]->answer) + write_contrast_colors(opt_output[o_intensity]->answer); + if (opt_output[o_exposition]->answer) + write_contrast_colors(opt_output[o_exposition]->answer); + if (opt_output[o_range]->answer) + write_contrast_colors(opt_output[o_range]->answer); } - G_percent(row, nrows, 2); /* end main loop */ - - /* finish and close */ - free_map(elevation.elev, row_buffer_size + 1); - for (i = o_forms; i < o_size; ++i) - if (opt_output[i]->answer) { - G_free(rasters[i].buffer); - Rast_close(rasters[i].fd); - Rast_short_history(opt_output[i]->answer, "raster", &history); - Rast_command_history(&history); - Rast_write_history(opt_output[i]->answer, &history); - } - - if (opt_output[o_forms]->answer) - write_form_cat_colors(opt_output[o_forms]->answer); - if (opt_output[o_intensity]->answer) - write_contrast_colors(opt_output[o_intensity]->answer); - if (opt_output[o_exposition]->answer) - write_contrast_colors(opt_output[o_exposition]->answer); - if (opt_output[o_range]->answer) - write_contrast_colors(opt_output[o_range]->answer); + G_free(strip_block); + G_free(strip_ptr); G_done_msg(" "); if (oneoff) { @@ -652,3 +777,50 @@ int main(int argc, char **argv) exit(EXIT_SUCCESS); } } + +static void compute_forms(FCELL **rows, int cur_row, int row, int col, + const struct geomorphon_config *cfg, + PATTERN *patterns, struct geomorphon_result *res) +{ + double eff_search = cfg->search_dist; + double eff_skip = cfg->skip_dist; + double eff_flat = cfg->flat_dist; + PATTERN *pattern; + int pattern_size; + FORMS cur_form, orig_form; + + pattern_size = calc_pattern(&patterns[0], row, cur_row, col, cfg->oneoff, + eff_search, eff_flat, rows); + pattern = &patterns[0]; + cur_form = orig_form = + determine_form(pattern->num_negatives, pattern->num_positives); + + if (cfg->extended && eff_search > 10 * cfg->max_resolution) { + if ((cur_form == SH || cur_form == FS || cur_form == PK || + cur_form == RI)) { + FORMS small_form; + + eff_search = (cfg->search_dist / 2. < 4 * cfg->max_resolution) + ? 4 * cfg->max_resolution + : cfg->search_dist / 2.; + eff_skip = 0; + eff_flat = 0; + pattern_size = calc_pattern(&patterns[1], row, cur_row, col, 0, + eff_search, eff_flat, rows); + pattern = &patterns[1]; + small_form = + determine_form(pattern->num_negatives, pattern->num_positives); + if (cur_form == SH || cur_form == FS) + cur_form = (small_form == FL) ? FL : cur_form; + if (cur_form == PK || cur_form == RI) + cur_form = small_form; + } + } + res->pattern = pattern; + res->pattern_size = pattern_size; + res->cur_form = cur_form; + res->orig_form = orig_form; + res->eff_search = eff_search; + res->eff_skip = eff_skip; + res->eff_flat = eff_flat; +} diff --git a/raster/r.geomorphon/memory.c b/raster/r.geomorphon/memory.c index b615772f81b..64aab12271a 100644 --- a/raster/r.geomorphon/memory.c +++ b/raster/r.geomorphon/memory.c @@ -35,11 +35,8 @@ static int get_cell(int, float *, void *, RASTER_MAP_TYPE); int open_map(MAPS *rast) { - - int row, col; char *mapset; struct Cell_head cellhd; - void *tmp_buf; mapset = (char *)G_find_raster2(rast->elevname, ""); @@ -57,17 +54,6 @@ int open_map(MAPS *rast) "Run g.region raster=%s to set proper resolution"), rast->elevname, rast->elevname); - tmp_buf = Rast_allocate_buf(rast->raster_type); - rast->elev = (FCELL **)G_malloc((row_buffer_size + 1) * sizeof(FCELL *)); - - for (row = 0; row < row_buffer_size + 1; ++row) { - rast->elev[row] = Rast_allocate_buf(FCELL_TYPE); - Rast_get_row(rast->fd, tmp_buf, row, rast->raster_type); - for (col = 0; col < ncols; ++col) - get_cell(col, rast->elev[row], tmp_buf, rast->raster_type); - } /* end elev */ - - G_free(tmp_buf); return 0; } @@ -102,38 +88,17 @@ static int get_cell(int col, float *buf_row, void *buf, return 0; } -int shift_buffers(int row) +/* Load count rows from abs_first into rows[], converted to FCELL. */ +int load_strip(int fd, RASTER_MAP_TYPE rtype, void *tmp_buf, FCELL **rows, + int abs_first, int count) { - int i; - int col; - void *tmp_buf; - FCELL *tmp_elev_buf; - - tmp_buf = Rast_allocate_buf(elevation.raster_type); - tmp_elev_buf = elevation.elev[0]; - - for (i = 1; i < row_buffer_size + 1; ++i) - elevation.elev[i - 1] = elevation.elev[i]; - - elevation.elev[row_buffer_size] = tmp_elev_buf; - Rast_get_row(elevation.fd, tmp_buf, row + row_radius_size + 1, - elevation.raster_type); - - for (col = 0; col < ncols; ++col) - get_cell(col, elevation.elev[row_buffer_size], tmp_buf, - elevation.raster_type); + int k, col; - G_free(tmp_buf); - return 0; -} - -int free_map(FCELL **map, int n) -{ - int i; - - for (i = 0; i < n; ++i) - G_free(map[i]); - G_free(map); + for (k = 0; k < count; ++k) { + Rast_get_row(fd, tmp_buf, abs_first + k, rtype); + for (col = 0; col < ncols; ++col) + get_cell(col, rows[k], tmp_buf, rtype); + } return 0; } diff --git a/raster/r.geomorphon/pattern.c b/raster/r.geomorphon/pattern.c index 5615941d889..1447bf5e55b 100644 --- a/raster/r.geomorphon/pattern.c +++ b/raster/r.geomorphon/pattern.c @@ -68,7 +68,8 @@ static int compare_multi(const double nadir_angle, const double zenith_angle, } int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, - const int oneoff) + const int oneoff, double search_distance, double flat_distance, + FCELL **rows) { /* calculate parameters of geomorphons and store it in the struct pattern */ int i, j, pattern_size = 0; @@ -83,7 +84,7 @@ int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, /* use distance calculation */ cur_northing = Rast_row_to_northing(row + 0.5, &window); cur_easting = Rast_col_to_easting(col + 0.5, &window); - center_height = elevation.elev[cur_row][col]; + center_height = rows[cur_row][col]; pattern->num_positives = 0; pattern->num_negatives = 0; pattern->positives = 0; @@ -108,8 +109,7 @@ int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, cur_row + j * nextr[i] > row_buffer_size - 1 || col + j * nextc[i] < 0 || col + j * nextc[i] > ncols - 1) continue; /* border: current cell is on the end of DEM */ - if (Rast_is_f_null_value( - &elevation.elev[cur_row + nextr[i]][col + nextc[i]])) + if (Rast_is_f_null_value(&rows[cur_row + nextr[i]][col + nextc[i]])) continue; /* border: next value is null, line-of-sight does not exists */ pattern_size++; /* line-of-sight exists, continue calculate visibility @@ -134,9 +134,8 @@ int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, col + j * nextc[i] < 0 || col + j * nextc[i] > ncols - 1) break; /* reached end of DEM (cols) or buffer (rows) */ - height = - elevation.elev[cur_row + j * nextr[i]][col + j * nextc[i]] - - center_height; + height = rows[cur_row + j * nextr[i]][col + j * nextc[i]] - + center_height; angle = atan2(height, cur_distance); if (angle > zenith_angle) { diff --git a/raster/r.geomorphon/r.geomorphon.md b/raster/r.geomorphon/r.geomorphon.md index 9deeac0fc99..8a9aa717bf1 100644 --- a/raster/r.geomorphon/r.geomorphon.md +++ b/raster/r.geomorphon/r.geomorphon.md @@ -195,6 +195,27 @@ elevation difference of at least 136 m is required to be noticed as non-flat. Flatness distance threshold may be helpful to avoid this problem. +### Performance + +To enable parallel processing, the user can specify the number of threads +with the **nprocs** parameter (default 1). The **memory** parameter +(default 300, in MB) sets how much memory the computation may use for its +buffers. The region is processed in horizontal bands whose height is +derived from this value, so a larger setting gives taller bands and fewer +passes over the input. Very small values still produce the same result but +re-read the search window halo more often. To take advantage of the +parallelization, GRASS needs to be compiled with OpenMP enabled. + +On an 8 core Apple M3, the forms output scales by roughly 1.9x on 2 +threads, 3.5x on 4 threads and 5.3x on 8 threads on maps from 10 to 100 +million cells. These numbers are medians from cooldown-controlled runs on +that machine; the benchmark script in the module source directory reports +averages instead. + +Actual memory use can slightly exceed the **memory** value because each +thread keeps a small per thread scratch row and the shared input strip +holds one extra row beyond the band and its halo. + ## EXAMPLES ### Geomorphon calculation: extraction of terrestrial landforms diff --git a/raster/r.geomorphon/tests/r_geomorphon_test.py b/raster/r.geomorphon/tests/r_geomorphon_test.py index 38d7e69546c..c8ce4b4c705 100644 --- a/raster/r.geomorphon/tests/r_geomorphon_test.py +++ b/raster/r.geomorphon/tests/r_geomorphon_test.py @@ -65,6 +65,22 @@ ("forms_skip1", "forms", {"skip": 1}), ] +# All twelve raster outputs, produced in one call by the parallel tests. +OUTPUTS = [ + "forms", + "ternary", + "positive", + "negative", + "intensity", + "exposition", + "range", + "variance", + "elongation", + "azimuth", + "extend", + "width", +] + # Landform categories from flat (1) to pit (10). ALL_CLASSES = set(range(1, 11)) @@ -163,3 +179,103 @@ def test_profile_json(fixed_region): assert inter["num_positives"] == 3 assert inter["num_negatives"] == 5 assert inter["pattern_size"] == 8 + + +def _assert_pairs_identical(tools, pairs): + """Check that the two rasters in each pair match everywhere, both in value and in which cells are null.""" + value_expr = " + ".join(f"abs(double({a}) - double({b}))" for a, b in pairs) + null_expr = " + ".join(f"(isnull({a}) != isnull({b}))" for a, b in pairs) + tools.r_mapcalc(expression=f"_vdiff = {value_expr}", overwrite=True) + tools.r_mapcalc(expression=f"_ndiff = {null_expr}", overwrite=True) + # A nonzero maximum in either combined map means the outputs diverged. + assert tools.r_univar(map="_ndiff", format="json")["max"] == 0 + assert tools.r_univar(map="_vdiff", format="json")["max"] == 0 + + +class TestParallelIdentity: + """Threaded output must equal serial output exactly.""" + + def _run_all(self, tools, suffix, nprocs, extra=None): + # One call produces all twelve outputs at the given thread count. + call = { + "elevation": "dem", + "search": SEARCH, + "nprocs": nprocs, + "overwrite": True, + } + for out in OUTPUTS: + call[out] = f"{out}_{suffix}" + if extra: + call.update(extra) + tools.r_geomorphon(**call) + + def test_threads_1_vs_4(self, fixed_region): + """All twelve outputs match between nprocs=1 and nprocs=4.""" + tools = Tools(session=fixed_region) + self._run_all(tools, "s", 1) + self._run_all(tools, "p", 4) + _assert_pairs_identical(tools, [(f"{o}_s", f"{o}_p") for o in OUTPUTS]) + + def test_band_seam_identity(self, fixed_region): + """memory=0 forces multiple bands and the output still matches serial.""" + tools = Tools(session=fixed_region) + self._run_all(tools, "s", 1) + self._run_all(tools, "b", 4, extra={"memory": 0}) + _assert_pairs_identical(tools, [(f"{o}_s", f"{o}_b") for o in OUTPUTS]) + + def test_min_legal_height(self, geomorphon_session): + """A region of exactly 2*search+2 rows at nprocs=4 matches serial.""" + tools = Tools(session=geomorphon_session) + tools.g_region(s=0, n=12, w=0, e=40, res=1) + tools.r_geomorphon( + elevation="dem", forms="forms_s12", search=SEARCH, nprocs=1, overwrite=True + ) + tools.r_geomorphon( + elevation="dem", forms="forms_p12", search=SEARCH, nprocs=4, overwrite=True + ) + _assert_pairs_identical(tools, [("forms_s12", "forms_p12")]) + + def test_mask_identity(self, fixed_region): + """With a mask active, threaded output matches serial.""" + tools = Tools(session=fixed_region) + tools.r_mapcalc(expression="msk = if(col() < 20, 1, null())", overwrite=True) + tools.r_mask(raster="msk") + try: + tools.r_geomorphon( + elevation="dem", + forms="forms_sm", + search=SEARCH, + nprocs=1, + overwrite=True, + ) + tools.r_geomorphon( + elevation="dem", + forms="forms_pm", + search=SEARCH, + nprocs=4, + overwrite=True, + ) + _assert_pairs_identical(tools, [("forms_sm", "forms_pm")]) + finally: + tools.r_mask(flags="r") + + def test_extended_identity(self, fixed_region): + """The extended correction path matches between serial and threaded.""" + tools = Tools(session=fixed_region) + tools.r_geomorphon( + elevation="dem", + forms="forms_se", + search=12, + flags="e", + nprocs=1, + overwrite=True, + ) + tools.r_geomorphon( + elevation="dem", + forms="forms_pe", + search=12, + flags="e", + nprocs=4, + overwrite=True, + ) + _assert_pairs_identical(tools, [("forms_se", "forms_pe")])