diff --git a/README.org b/README.org index 32ed854..dcac33b 100644 --- a/README.org +++ b/README.org @@ -127,9 +127,46 @@ bounds, run the =./standalone= tool. Example: [[file:example-standalone.png]] This can either run a static GLUT application, or it can render to a =.png= -image on disk and/or a binary range image. Run with =--help= for details. Note -that the azimuth extents are currently specified differently than they are in -the interactive tool. +image on disk (or a =.pdf=/=.svg= image annotated with named peaks; see +below) and/or a binary range image. Note that the azimuth extents are +currently specified differently than they are in the interactive tool. + +*** Options +Run =./standalone --help= for the full, up-to-date list of options and a +detailed explanation of how each one affects the render. The most important +ones, grouped by what they do: + +- *Field of view*: =--znear=/=--zfar= bound the rendered distance, in + meters; raising =--znear= is a good way to skip an uninteresting nearby + plain. The vertical field of view isn't set directly: it falls out of + =--width=, =--height= and =AZ_RADIUS_DEG=, since the same angular + resolution is used in both directions. A real mountain panorama usually + only occupies a few degrees of elevation, so an explicit, small + =--height= (rather than the default, which assumes a much wider vertical + span) both crops away the resulting excess sky and increases the + effective resolution of the relief that matters -- the single biggest + lever for a sharp-looking render. +- *Earth curvature and atmospheric refraction*: =--curvature= corrects the + apparent elevation of distant terrain for the curvature of the Earth, + partly compensated by atmospheric refraction (=--refraction-k=, default + 0.13). Off by default (flat tangent-plane rendering); matters at 100+km. +- *Color and ridge outlines*: by default the render is color-coded by + range (atmospheric perspective: dark neutral gray near, fading to a pale + blue-gray far away, closer to the white background), with the crest of + each visible ridge/mountain-range layer outlined, the way + classic drawn panoramas (e.g. [[https://www.udeuschle.de/panoramas/makepanoramas_en.htm][udeuschle.de]]) do it. See + =--no-ridge-lines=, =--ridge-line-threshold= and =--ridge-line-gray= to + disable or tune this. +- *Peak labels*: with a =.pdf=/=.svg= output filename, named peaks are + labelled on the render (see [[https://github.com/dkogan/horizonator/blob/master/query-peaks-from-osm.py][query-peaks-from-osm.py]] to generate the list + for your own area, and edit the =#include= near the bottom of + [[https://github.com/dkogan/horizonator/blob/master/standalone.c][standalone.c]] to point at it). +- *Viewer position*: =--viewer-height= adds a number of meters to the + viewer elevation sampled from the DEM (e.g. the height of an apartment + floor above street level). +- *Performance*: by default the mesh only covers the azimuth wedge that's + actually being rendered, not the full circle of loaded DEM data; see + =--no-restrict-mesh-azimuth= and the DEM resolution section below. ** C API The tool can be invoked from C. The [[https://github.com/dkogan/horizonator/blob/master/horizonator.h][header comments]] and its usages in the @@ -155,10 +192,13 @@ things works well. The view straight ahead (elevation = 0) is at the center of the render. -This tool operates in the tangent plane to the viewer, so it assumes that -locally, the Earth is flat. This produces small inaccuracies, but unless we care -about small pixel-level errors, this is a good approximation. I will eventually -fix this. +By default this tool operates in the tangent plane to the viewer, so it +assumes that locally, the Earth is flat. This is a good approximation at +short range, but at 100+km (e.g. distant high peaks) the apparent elevation +angle can be off by upwards of a kilometer of apparent height. Pass +=--curvature= (on the =standalone= tool; =horizonator_set_curvature()= in +the C/Python APIs) to correct for Earth curvature and atmospheric +refraction instead. * DEM resolution By default, 3" DEMs are used. These are the low-res SRTM data, which is @@ -172,18 +212,30 @@ objects. Until that is implemented, the 9x increase in triangles present in the SRTM1 data could become a problem. Support for 1" data /is/ in place, and can be selected with the =SRTM1= option in all the APIs and commandline tools. +The =standalone= tool mitigates this somewhat: since each invocation renders a +single, fixed azimuth wedge, it only meshes that wedge (plus a small margin) +instead of the full circle of loaded DEM data, cutting the triangle count +without changing the output (see =--no-restrict-mesh-azimuth= to disable +this). This doesn't help the interactive =horizonator= tool, which allows +panning to any azimuth after the data is loaded, and so still meshes the +full circle. + * Nice-to-have improvements In no particular order: - Texturing with aerial imagery - Being more efficient about data loading: the DEM and texture resolution needs - to be high close-in, but can be dramatically lower further out. + to be high close-in, but can be dramatically lower further out. The + =standalone= tool only meshes the azimuth wedge it renders (see DEM + resolution above), but a true level-of-detail scheme (full resolution + close in, coarser far away) isn't implemented yet - Higher-res DEMs are available (1sec SRTM instead of 3sec). It would be nice to use them, /if/ we can do so efficiently - Nicer handling of the mesh immediately near the viewer. - Intelligently loading faraway data. Currently we load data a constant number of cells away from the viewer -- Peak-labelling the render +- Peak-labelling the render: done for the =standalone= tool's =.pdf=/=.svg= + output (see [[https://github.com/dkogan/horizonator/blob/master/query-peaks-from-osm.py][query-peaks-from-osm.py]]); not yet wired into the interactive tool - More UI stuff - text showing the current lat, lon, az bounds - text inputs to change the current lat, lon, az bounds diff --git a/annotator.c b/annotator.c index 1ef918c..9df26b1 100644 --- a/annotator.c +++ b/annotator.c @@ -16,22 +16,55 @@ #include "horizonator.h" -#define MAX_MARKER_DIST 100000.0 #define MIN_MARKER_DIST 500.0 #define FUZZ_RANGE 500. -#define FUZZ_PIXEL_Y 6 + +// The vertical pixel search radius used to find a POI's true rendered +// position isn't a fixed pixel count: the real positional uncertainty +// (OSM coordinate precision, DEM sampling) is an ANGULAR one, a small +// fraction of a degree. At low resolution a handful of pixels happens to +// cover that; at high resolution the same few pixels cover only a tiny +// sliver of a degree, and otherwise-valid matches silently fail. So a +// fixed angular tolerance is converted to a pixel count from the actual +// image resolution (see fuzz_pixel_y() below) instead +#define ANGULAR_FUZZ_DEG 0.05 #define LABEL_CROSSHAIR_R 3 #define TEXT_MARGIN 2 +// Label text in black, leader line in a soft "almond green" -- readable on +// the white background, and the line doesn't compete visually with the text +#define LABEL_TEXT_R 0.0 +#define LABEL_TEXT_G 0.0 +#define LABEL_TEXT_B 0.0 +#define LABEL_LINE_R 0.576 +#define LABEL_LINE_G 0.773 +#define LABEL_LINE_B 0.447 + +// Label text is drawn rotated by this angle: first letter at the bottom, +// last letter at the top (a vertical column of text, reading bottom-up, +// as if tilting your head to the left). cairo_rotate() rotates towards +// the (downward-pointing) +y axis for positive angles, so a negative +// angle here goes up, as wanted +#define LABEL_ROTATION_RAD (-M_PI/2.0) + +// Every kept label's last (topmost) character sits this many pixels below +// the top edge of the canvas (not the rendered/mesh area -- the canvas) +#define LABEL_TOP_MARGIN_PX 30.0 + +// Two labels whose horizontal (screen-x) positions are closer than this +// are considered to conflict, and only the higher-elevation one of the +// group is kept. Since the label text is vertical, its own on-screen +// footprint is about one font-height wide; this multiplies that by a +// safety margin +#define LABEL_CONFLICT_GAP_FACTOR 1.5 + static const double POINTS_PER_INCH = 72.; static const double PIXELS_PER_INCH = 300.; static const double CAIRO_SCALE = POINTS_PER_INCH / PIXELS_PER_INCH; -static int font_height = 20; - static double string_width(cairo_t *cr, const char* s) @@ -47,19 +80,13 @@ double string_width(cairo_t *cr, -typedef struct -{ - float x,y; -} xy_t; - -// compares two POIs by their draw_x. Sorts disabled POIs to the end -static int compar_poi_x( const void* _idx0, const void* _idx1, void* cookie ) +// compares two visible_poi_t by their screen-x position +static int compar_visible_poi_x( const void* _a, const void* _b ) { - const xy_t* xy = (const xy_t*)cookie; - const int* idx0 = (const int*)_idx0; - const int* idx1 = (const int*)_idx1; + const visible_poi_t* a = (const visible_poi_t*)_a; + const visible_poi_t* b = (const visible_poi_t*)_b; - if( xy[ *idx0 ].x < xy[ *idx1 ].x ) + if( a->x < b->x ) return -1; else return 1; @@ -67,23 +94,42 @@ static int compar_poi_x( const void* _idx0, const void* _idx1, void* cookie ) static void draw_label( cairo_t* cr, + // the peak's own screen position double x, double y, - // top of the label - double y_label, double lat, double lon, const char* name ) { + // The text is vertical (LABEL_ROTATION_RAD), growing from its first + // (bottom) character up to its last (top) character. We want the LAST + // character LABEL_TOP_MARGIN_PX below the top of the canvas, so the + // anchor (first character, where the leader line ends) is that margin + // plus the full string length below the canvas top + const double string_w = string_width(cr, name); + const double text_start_y = LABEL_TOP_MARGIN_PX + string_w; + + cairo_set_source_rgb(cr, LABEL_LINE_R, LABEL_LINE_G, LABEL_LINE_B); + cairo_move_to(cr, x-LABEL_CROSSHAIR_R, y); cairo_rel_line_to(cr, 2*LABEL_CROSSHAIR_R, 0); cairo_move_to(cr, x, y+LABEL_CROSSHAIR_R); - cairo_line_to(cr, x, y_label); + cairo_line_to(cr, x, text_start_y); cairo_stroke(cr); + // The label text, rotated LABEL_ROTATION_RAD: anchored exactly where the + // leader line ends, and drawn along the rotated axis from there, so the + // anchor point doesn't move as the rotation changes. cairo's "current + // point" is a user-space coordinate, so to rotate around a specific + // device-space point we translate there FIRST, rotate, then draw at the + // new (local) origin + cairo_set_source_rgb(cr, LABEL_TEXT_R, LABEL_TEXT_G, LABEL_TEXT_B); + + cairo_save(cr); + cairo_translate(cr, x, text_start_y); + cairo_rotate(cr, LABEL_ROTATION_RAD); + cairo_move_to(cr, 0, 0); - // cairo wants the bottom of the label - cairo_move_to(cr, x, y_label + font_height); char url[256]; bool url_valid = (snprintf(url, sizeof(url), @@ -93,6 +139,8 @@ void draw_label( cairo_t* cr, if(url_valid) cairo_tag_begin (cr, CAIRO_TAG_LINK, url); cairo_show_text(cr, name); if(url_valid) cairo_tag_end (cr, CAIRO_TAG_LINK); + + cairo_restore(cr); } @@ -139,6 +187,112 @@ bool RGB32_from_BGR24(// output return result; } +int find_visible_pois(// output + visible_poi_t* visible, + + // input + const float* range_image, + const int width, + const int height, + const int cut_off_bottom_px, + + const poi_t* pois, + const int Npois, + const double lat, + const double lon, + const double az_deg0, + const double az_deg1, + const double ele_m, + + const bool curvature_enabled, + const double refraction_k, + + const double max_marker_dist_m) +{ + const int height_out = height - cut_off_bottom_px; + const double cos_lat = cos(lat * M_PI/180.); + + // Same angular resolution horizontally and vertically (see README), so + // this is the degrees/pixel in both directions + const double deg_per_pixel = (az_deg1-az_deg0) / (double)width; + const int fuzz_pixel_y = (int)ceil(ANGULAR_FUZZ_DEG / deg_per_pixel); + + int Nvisible = 0; + + for(int i=0; i max_marker_dist_m ) + // too close or too far to label + continue; + + // I'm finished with the projection. I now unproject to look for + // occlusions + + // The rendered peaks usually don't end up exactly where the POI list + // says they should be. I scan the range map vertically to find the true + // peak (or to decide that it's occluded) + int fuzz_nearest = 0; // initializing to pacify compiler + double err_nearest = DBL_MAX; + + for( int fuzz = -fuzz_pixel_y; fuzz < fuzz_pixel_y; fuzz++ ) + { + if(crosshair_y + (double)fuzz < 0) + continue; + if( crosshair_y + (double)fuzz >= height_out ) + break; + + // As I move down the image the range will get closer and closer. I + // pick the highest value that's closest + const float range = + range_image[width*( (int)round(crosshair_y) + fuzz) + + (int)round(crosshair_x)]; + + if(range <= 0.0f) + // no render data here + continue; + + double err = fabs(range_have - range); + if( err < err_nearest ) + { + err_nearest = err; + fuzz_nearest = fuzz; + } + else + // it can only get worse from here, so give up + break; + } + + if( err_nearest < FUZZ_RANGE ) + { + visible[Nvisible].poi_index = i; + visible[Nvisible].x = crosshair_x; + visible[Nvisible].y = crosshair_y + (double)fuzz_nearest; + visible[Nvisible].range = range_have; + Nvisible++; + } + } + + return Nvisible; +} + bool annotate(// input const char* out_filename, // assumed to be stored densely. @@ -154,17 +308,33 @@ bool annotate(// input const double lon, const double az_deg0, const double az_deg1, - const double ele_m) + const double ele_m, + + const bool curvature_enabled, + const double refraction_k, + + // POIs farther than this are never labelled. Pass the same + // zfar used for the render: there's no point labelling + // something farther than what was actually rendered + const double max_marker_dist_m, + + // Label text height, in real typographic points (1/72in), + // as it'll appear in the output PDF/SVG + const double label_font_size_pt) { bool result = false; - const int height_out = height - cut_off_bottom_px; + // label_font_size_pt is in real typographic points, i.e. as it appears + // in the final PDF/SVG page (which is CAIRO_SCALE units per image + // pixel). font_height is the corresponding size in the image's own + // pixel-like coordinate system, which is what cairo_set_font_size() + // wants here, since we're drawing under a cairo_scale(CAIRO_SCALE) + const double font_height = label_font_size_pt / CAIRO_SCALE; - // For sorting, further down - int poi_indices[Npois]; - int Npoi_indices = 0; + const int height_out = height - cut_off_bottom_px; - xy_t labels_xy[Npois]; + visible_poi_t visible[Npois]; + int Nvisible = 0; uint8_t* image_rgb32 = NULL; cairo_surface_t* surface = NULL; @@ -273,120 +443,68 @@ bool annotate(// input cairo_paint(cr); cairo_set_font_size(cr, font_height - TEXT_MARGIN); - cairo_set_source_rgb(cr, 1.0, 1.0, 0.0); - ////// Pick and render the annotations - for(int i=0; i MAX_MARKER_DIST ) - // too close or too far to label - continue; - - // crosshair_y will be checked below in the fuzz loop + qsort( visible, Nvisible, sizeof(visible[0]), &compar_visible_poi_x ); + const double conflict_gap = font_height * LABEL_CONFLICT_GAP_FACTOR; - // I'm finished with the projection. I now unproject to look for - // occlusions - - // The rendered peaks usually don't end up exactly where the POI list - // says they should be. I scan the range map vertically to find the true - // peak (or to decide that it's occluded) - int fuzz_nearest = 0; // initializing to pacify compiler - double err_nearest = DBL_MAX; + bool keep[Nvisible ? Nvisible : 1]; - for( int fuzz = -FUZZ_PIXEL_Y; fuzz < FUZZ_PIXEL_Y; fuzz++ ) - { - if(crosshair_y + (double)fuzz < 0) - continue; - if( crosshair_y + (double)fuzz >= height_out ) - break; - - // As I move down the image the range will get closer and closer. I - // pick the highest value that's closest - const float range = - range_image[width*( (int)round(crosshair_y) + fuzz) + - (int)round(crosshair_x)]; - - if(range <= 0.0f) - // no render data here - continue; - - double err = fabs(range_have - range); - if( err < err_nearest ) - { - err_nearest = err; - fuzz_nearest = fuzz; - } - else - // it can only get worse from here, so give up - break; - } - - if( err_nearest < FUZZ_RANGE ) - { - poi_indices[Npoi_indices++] = i; - labels_xy[i].x = crosshair_x; - labels_xy[i].y = crosshair_y + (float)fuzz_nearest; - } + int i = 0; + while( i < Nvisible ) + { + // Chain together consecutive (in x order) POIs whose gap to their + // neighbor is under the threshold: [i,j] is one conflicting group + int j = i; + while( j+1 < Nvisible && + visible[j+1].x - visible[j].x < conflict_gap ) + j++; + + int ibest = i; + for( int k=i+1; k<=j; k++ ) + if( pois[ visible[k].poi_index ].ele_m > pois[ visible[ibest].poi_index ].ele_m ) + ibest = k; + + for( int k=i; k<=j; k++ ) + keep[k] = (k == ibest); + + i = j+1; } - // Now that I have all the crosshair positions, compute the label positions. - - // start out by sorting the POIs by their crosshair_x - qsort_r( poi_indices, Npoi_indices, sizeof(poi_indices[0]), - &compar_poi_x, labels_xy ); - - // I now traverse the sorted list of POIs, keeping track of groups of POIs - // that overlap in the horizontal. After each overlapping group is complete, - // set up the labels of each group member to stagger the labels and avoid - // overlap - float overlapgroup_right = -1; // not in an overlapping group at first - float current_y = 0; // start on top - for( int i=0; ix; - float right = label_xy->x + string_width(cr,poi->name); - - if( left > overlapgroup_right || current_y + font_height >= height_out ) - { - // not overlapping, or the label is too low. Draw label on top. - current_y = 0; - overlapgroup_right = right; - } - else - { - // I overlap the previous. Thus draw the label a bit lower - if( overlapgroup_right < right ) - overlapgroup_right = right; - } + const poi_t* poi = &pois[ visible[i].poi_index ]; draw_label(cr, - label_xy->x, label_xy->y, current_y, + visible[i].x, visible[i].y, poi->lat, poi->lon, poi->name); - - current_y += font_height; } + } + + cairo_set_source_rgb(cr, LABEL_TEXT_R, LABEL_TEXT_G, LABEL_TEXT_B); const int bearing_annotation_spacing = 15; for(int az=180; az>-180; az -= bearing_annotation_spacing) diff --git a/annotator.h b/annotator.h index 73a6293..014e680 100644 --- a/annotator.h +++ b/annotator.h @@ -7,6 +7,46 @@ typedef struct float lat, lon, ele_m; } poi_t; +typedef struct +{ + int poi_index; // index into the pois[] array that was passed in + double x,y; // screen position in the render + double range; // distance from the viewer, in meters +} visible_poi_t; + +// Figures out which of pois[] are actually visible in the given render (not +// occluded by closer terrain), same logic annotate() uses internally to +// place labels. Returns the count of visible POIs, with that many entries +// filled in in visible[] (which must have room for Npois entries) +int find_visible_pois(// output: room for Npois entries + visible_poi_t* visible, + + // input + const float* range_image, + const int width, + const int height, + const int cut_off_bottom_px, + + const poi_t* pois, + const int Npois, + const double lat, + const double lon, + const double az_deg0, + const double az_deg1, + const double ele_m, + + // Must match whatever was passed to horizonator_set_curvature() + // for this render, or POIs will be placed incorrectly and fail + // to be found in the range image (looking occluded) + const bool curvature_enabled, + const double refraction_k, + + // POIs farther than this are never returned. Pass the + // same zfar used for the render: there's no point + // considering something farther than what was + // actually rendered + const double max_marker_dist_m); + bool annotate(// input const char* out_filename, // must be .pdf or .svg // assumed to be stored densely. @@ -22,4 +62,19 @@ bool annotate(// input const double lon, const double az_deg0, const double az_deg1, - const double ele_m); + const double ele_m, + + // Must match whatever was passed to horizonator_set_curvature() + // for this render, or POIs will be placed incorrectly and fail + // to be found in the range image (looking occluded) + const bool curvature_enabled, + const double refraction_k, + + // POIs farther than this are never labelled. Pass the same + // zfar used for the render: there's no point labelling + // something farther than what was actually rendered + const double max_marker_dist_m, + + // Label text height, in real typographic points (1/72in), + // as it'll appear in the output PDF/SVG + const double label_font_size_pt); diff --git a/cluster-visible-peaks.py b/cluster-visible-peaks.py new file mode 100755 index 0000000..bb2bcb8 --- /dev/null +++ b/cluster-visible-peaks.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python3 +"""Group the peaks visible in a render into "massifs" by geographic proximity + +Reads the tab-separated output of + + ./standalone --list-visible-peaks OUT.txt ... + +(one line per visible peak: name, lat, lon, ele_m, range_m), drops peaks +with no letters in their name (OSM had no real name for these; the query +script falls back to printing their altitude instead, e.g. "1583.0"), and +groups the rest into "massifs" via simple distance-based clustering: any +two peaks closer than --cluster-km are considered part of the same group +(transitively, so a chain of nearby peaks can span a wider area than +--cluster-km end to end). + +There's no reliable "massif" tag in OpenStreetMap for this area, so these +groups have no real name -- they're just numbered, ordered by mean +distance from the viewer (nearest massif first). Rename/relabel them by +hand afterwards if wanted. + +Usage: + ./cluster-visible-peaks.py visible-peaks.txt [--cluster-km 8] +""" + +import argparse +import math +import sys + + +def haversine_km(lat1, lon1, lat2, lon2): + R = 6371.0 + p1, p2 = math.radians(lat1), math.radians(lat2) + dphi = math.radians(lat2 - lat1) + dlambda = math.radians(lon2 - lon1) + a = math.sin(dphi/2)**2 + math.cos(p1)*math.cos(p2)*math.sin(dlambda/2)**2 + return 2*R*math.asin(math.sqrt(a)) + + +def has_letter(s): + return any(c.isalpha() for c in s) + + +def read_peaks(filename): + peaks = [] + with open(filename, encoding='utf-8') as f: + for line in f: + line = line.rstrip('\n') + if not line: + continue + name, lat, lon, ele, rng = line.split('\t') + peaks.append({'name': name, + 'lat': float(lat), + 'lon': float(lon), + 'ele': float(ele), + 'range': float(rng)}) + return peaks + + +# Union-find: transitively chains together any peaks closer than +# cluster_km, same idea as the horizontal-conflict grouping in annotator.c +def cluster(peaks, cluster_km): + n = len(peaks) + parent = list(range(n)) + + def find(x): + while parent[x] != x: + parent[x] = parent[parent[x]] + x = parent[x] + return x + + def union(a, b): + ra, rb = find(a), find(b) + if ra != rb: + parent[ra] = rb + + for i in range(n): + for j in range(i+1, n): + if haversine_km(peaks[i]['lat'], peaks[i]['lon'], + peaks[j]['lat'], peaks[j]['lon']) < cluster_km: + union(i, j) + + groups = {} + for i in range(n): + groups.setdefault(find(i), []).append(i) + return groups + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('infile', + help="Output of standalone --list-visible-peaks") + parser.add_argument('--cluster-km', type=float, default=8.0, + help="Peaks closer than this (km) are grouped into the same massif (default 8)") + args = parser.parse_args() + + peaks = read_peaks(args.infile) + n_before = len(peaks) + peaks = [p for p in peaks if has_letter(p['name'])] + + print(f"# {n_before} sommets visibles, {len(peaks)} avec un vrai nom " + f"(altitude-only exclus)", file=sys.stderr) + + groups = cluster(peaks, args.cluster_km) + + def massif_mean_range(idxs): + return sum(peaks[i]['range'] for i in idxs) / len(idxs) + + ordered = sorted(groups.values(), key=massif_mean_range) + + for mi, idxs in enumerate(ordered): + idxs = sorted(idxs, key=lambda i: peaks[i]['range']) + mean_km = massif_mean_range(idxs) / 1000. + print(f"\n=== Massif {mi+1} ({len(idxs)} sommets, ~{mean_km:.0f} km) ===") + for i in idxs: + p = peaks[i] + print(f" {p['name']:<45s} ele={p['ele']:6.0f}m " + f"dist={p['range']/1000:5.1f}km " + f"({p['lat']:.5f}, {p['lon']:.5f})") + + +if __name__ == '__main__': + main() diff --git a/fragment.glsl b/fragment.glsl index 2882e5d..cb4fa30 100644 --- a/fragment.glsl +++ b/fragment.glsl @@ -3,21 +3,100 @@ #version 420 layout(location = 0) out vec4 frag_color; -in vec3 rgb_fragment; in vec2 tex_fragment; uniform sampler2D tex; uniform int NtilesX, NtilesY; +// shading_scale is 0.0 (disabled: legacy unshaded rendering) or 1.0 +// (enabled). sun_dir is a unit vector in the same (east,north,height) +// frame as normal_fragment, pointing FROM the terrain TOWARDS the sun +uniform float shading_scale; +uniform vec3 sun_dir; +in vec3 normal_fragment; + +// materials_scale is 0.0 (disabled: no material tint) or 1.0 (enabled). +// A first, purely procedural land-cover approximation -- no aerial +// imagery or land-cover data involved -- classifying each fragment by +// elevation (snow line) and slope steepness (bare rock), with forest +// below the tree line and alpine grass above it otherwise +uniform float materials_scale; +in float elevation_fragment; + +// 0 at znear_color, 1 at zfar_color (see vertex.glsl) +in float atmo_t_fragment; + +const float SNOW_LINE_M = 2800.0; // above this: snow, any slope +const float TREE_LINE_M = 1800.0; // below (and not snow/rock): forest +const float SLOPE_ROCK_NZ = 0.55; // normal.z under this: bare rock + +const vec3 COLOR_FOREST = vec3(0.35, 0.42, 0.30); +const vec3 COLOR_GRASS = vec3(0.55, 0.58, 0.38); +const vec3 COLOR_ROCK = vec3(0.50, 0.47, 0.43); +const vec3 COLOR_SNOW = vec3(0.95, 0.96, 0.98); + +// Used as the near (atmo_t==0) color when materials_scale==0: this is +// the plain-grayscale look (no material data), a dark neutral gray purely +// for visual weight/contrast -- not a claim about what's really there. +// When materials ARE on, the material's own true color is used as the +// near color instead (see main() below): a real land-cover color +// shouldn't ALSO be darkened by this legacy gray, or the two darkening +// effects compound (near-dark-gray times a dim material times slope +// shading was crushing near, shaded forest to near black) +const vec3 COLOR_NEAR_DEFAULT = vec3(0.30, 0.30, 0.30); + +// The far (atmo_t==1) color, regardless of materials: a pale blue-gray +// haze, the way atmospheric scattering tints distant relief in a real +// photo. This one genuinely is distance-dependent physics (haze), unlike +// COLOR_NEAR_DEFAULT above, so it applies the same whether or not +// materials are on +const vec3 COLOR_FAR = vec3(0.80, 0.84, 0.92); + +vec3 material_color(float elevation, float slope_nz) +{ + if(elevation > SNOW_LINE_M) + return COLOR_SNOW; + if(slope_nz < SLOPE_ROCK_NZ) + return COLOR_ROCK; + if(elevation > TREE_LINE_M) + return COLOR_GRASS; + return COLOR_FOREST; +} void main(void) { + // normal_fragment is linearly interpolated (by the rasterizer) from + // the 3 per-vertex normals of the enclosing triangle, so it isn't + // unit length in general; renormalize before using it + vec3 n = normalize(normal_fragment); + + float ndotl = max(dot(n, sun_dir), 0.0); + // Ambient floor of 0.5: a slope facing away from the sun is dimmer, + // not black + float light = mix(0.5, 1.0, ndotl); + // shading_scale==0 must reproduce the old, unshaded look exactly + float light_scaled = mix(1.0, light, shading_scale); + + // materials_scale==0 must reproduce the old, untinted look exactly: + // mix(COLOR_NEAR_DEFAULT, material_color(...), 0.0) == COLOR_NEAR_DEFAULT + vec3 near_color = mix(COLOR_NEAR_DEFAULT, + material_color(elevation_fragment, n.z), + materials_scale); + + // Atmospheric haze: blend towards the far color with distance. This + // is the ONLY distance-dependent darkening/tinting left (on top of + // whichever near_color was picked above) -- unlike before, materials + // no longer also get multiplied by a separate near/far gray gradient + vec3 rgb = mix(near_color, COLOR_FAR, atmo_t_fragment); + + vec3 rgb_shaded = rgb * light_scaled; + if(NtilesX == 0) - frag_color = vec4(rgb_fragment, 1.0); + frag_color = vec4(rgb_shaded, 1.0); else { vec4 texcolor = texture( tex, tex_fragment.xy); - vec4 shadingcolor = vec4(rgb_fragment, 0.0); + vec4 shadingcolor = vec4(rgb_shaded, 0.0); frag_color = 0.7*texcolor + 0.3*shadingcolor; } } diff --git a/geometry.glsl b/geometry.glsl index 17ee1e6..5da2cc2 100644 --- a/geometry.glsl +++ b/geometry.glsl @@ -5,11 +5,27 @@ layout (triangles) in; layout (triangle_strip, max_vertices=3) out; -in vec3 rgb[]; -out vec3 rgb_fragment; +// 0 at znear_color, 1 at zfar_color -- how far towards the atmospheric +// haze color this point should be blended (see COLOR_NEAR_DEFAULT in +// fragment.glsl for why that blend happens there, not here) +in float atmo_t[]; +out float atmo_t_fragment; in vec2 tex[]; out vec2 tex_fragment; +// Per-vertex normal (world-space east/north/height frame), from +// vertex.glsl. Just passed through here, per vertex, so the fragment +// shader gets a value smoothly interpolated by the rasterizer across each +// triangle -- and so continuous across the edge between two triangles +// that share a vertex, unlike a flat per-triangle normal +in vec3 normal[]; +out vec3 normal_fragment; + +// Raw DEM elevation, for the procedural material classification in +// fragment.glsl +in float elevation_m[]; +out float elevation_fragment; + void main() { // The azimuth is gl_Position.x. Any triangles on the seam (some vertices @@ -28,9 +44,11 @@ void main() for(int i=0; i<3; i++) { - rgb_fragment = rgb[i]; - tex_fragment = tex[i]; - gl_Position = gl_in[i].gl_Position; + atmo_t_fragment = atmo_t[i]; + tex_fragment = tex[i]; + normal_fragment = normal[i]; + elevation_fragment = elevation_m[i]; + gl_Position = gl_in[i].gl_Position; EmitVertex(); } EndPrimitive(); diff --git a/horizonator-lib.c b/horizonator-lib.c index 552bc14..002bb21 100644 --- a/horizonator-lib.c +++ b/horizonator-lib.c @@ -29,6 +29,16 @@ #define OSM_TILE_TEXTURE_NAME_DEFAULT "mapnik" #define OSM_TILE_TEXTURE_URL_FMT_DEFAULT "https://tile.openstreetmap.org/%d/%d/%d.png" +// Used by restrict_mesh_azimuth (see horizonator_init()). Cells beyond this +// azimuth margin outside [mesh_az_deg0,mesh_az_deg1] are skipped when +// building the mesh +#define MESH_AZIMUTH_MARGIN_DEG 5.0f +// Cells within this many DEM cells of the viewer are always meshed, +// regardless of azimuth: right next to the viewer, azimuth changes very +// quickly from one cell to the next, so an azimuth-only test is unreliable +// there. This inner disk is cheap to mesh in full regardless +#define MESH_INNER_RADIUS_CELLS 100 + #define assert_opengl() \ do { \ int error = glGetError(); \ @@ -75,6 +85,9 @@ bool horizonator_init( // output int render_radius_cells, // This should be given >0 float render_radius_m, // or this, but not both + bool restrict_mesh_azimuth, + float mesh_az_deg0, float mesh_az_deg1, + bool use_glut, bool render_texture, bool SRTM1, @@ -185,7 +198,7 @@ bool horizonator_init( // output glEnable(GL_DEPTH_TEST); glEnable(GL_CULL_FACE); - glClearColor(0, 0, 1, 0); + glClearColor(1, 1, 1, 0); if( !horizonator_dem_init( &ctx->dems, viewer_lat, viewer_lon, @@ -439,7 +452,31 @@ bool horizonator_init( // output GLshort* vertices = glMapBuffer(GL_ARRAY_BUFFER, GL_WRITE_ONLY); #endif + // Per-vertex normal, for smooth (Gouraud-style) slope shading: the + // rasterizer interpolates this across each triangle, so shading is + // continuous across triangle edges (no visible facets), unlike a + // flat per-triangle normal. Estimated by finite differences on the + // 4 DEM neighbors -- needs real (not integer-quantized) precision, + // so this is its own float VBO rather than packed into the + // position one above + GLuint normalBufID; + glGenBuffers(1, &normalBufID); + glBindBuffer(GL_ARRAY_BUFFER, normalBufID); + glEnableVertexAttribArray(1); + glBufferData(GL_ARRAY_BUFFER, Nvertices*3*sizeof(GLfloat), NULL, GL_STATIC_DRAW); + glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, NULL); + GLfloat* normals = glMapBuffer(GL_ARRAY_BUFFER, GL_WRITE_ONLY); + + // meters/cell, East-West and North-South. Same formula as the + // (disabled) CPU-side paths below, factored out since the normal + // computation needs it on every vertex + const float Rearth = 6371000.0f; + const float cos_viewer_lat_here = cosf( M_PI / 180.0f * viewer_lat ); + const float cellsize_ns = Rearth * (float)(M_PI/180.0) / (float)ctx->dems.cells_per_deg; + const float cellsize_ew = cellsize_ns * cos_viewer_lat_here; + int vertex_buf_idx = 0; + int normal_buf_idx = 0; for( int j=0; j<2*render_radius_cells; j++ ) { @@ -485,14 +522,112 @@ bool horizonator_init( // output vertices[vertex_buf_idx++] = j; vertices[vertex_buf_idx++] = z; #endif + + // Finite-difference normal from the 4 DEM neighbors + // (clamped at the edges of the loaded grid, where a + // neighbor is missing: falls back to a one-sided + // difference there instead of a centered one) + int i_prev = i>0 ? i-1 : i; + int i_next = i<2*render_radius_cells-1 ? i+1 : i; + int j_prev = j>0 ? j-1 : j; + int j_next = j<2*render_radius_cells-1 ? j+1 : j; + + int32_t z_i_prev = horizonator_dem_sample(&ctx->dems, i_prev, j); + int32_t z_i_next = horizonator_dem_sample(&ctx->dems, i_next, j); + int32_t z_j_prev = horizonator_dem_sample(&ctx->dems, i, j_prev); + int32_t z_j_next = horizonator_dem_sample(&ctx->dems, i, j_next); + + // Tangent vectors along the East and North grid directions + // (in meters), and the cross product of the two (East x + // North = Up, in a right-handed ENU frame -- so this always + // comes out pointing "up", no sign ambiguity to resolve, + // unlike a per-triangle normal from arbitrarily-wound + // vertices) + float tangent_east_x = (float)(i_next - i_prev) * cellsize_ew; + float tangent_east_z = (float)(z_i_next - z_i_prev); + float tangent_north_y = (float)(j_next - j_prev) * cellsize_ns; + float tangent_north_z = (float)(z_j_next - z_j_prev); + + float nx = -tangent_east_z * tangent_north_y; + float ny = -tangent_east_x * tangent_north_z; + float nz = tangent_east_x * tangent_north_y; + + float ninv = 1.0f / sqrtf(nx*nx + ny*ny + nz*nz); + normals[normal_buf_idx++] = nx*ninv; + normals[normal_buf_idx++] = ny*ninv; + normals[normal_buf_idx++] = nz*ninv; } } int res = glUnmapBuffer(GL_ARRAY_BUFFER); assert( res == GL_TRUE ); + assert( normal_buf_idx == Nvertices*3 ); + + glBindBuffer(GL_ARRAY_BUFFER, vertexBufID); + res = glUnmapBuffer(GL_ARRAY_BUFFER); + assert( res == GL_TRUE ); assert( vertex_buf_idx == Nvertices*3 ); } + // Optional azimuth restriction: precompute, for each DEM cell, whether + // it lies within [mesh_az_deg0,mesh_az_deg1] (plus a margin), or close + // enough to the viewer to always be included. NULL means "no + // restriction: mesh the whole loaded circle", to keep the existing + // behavior for callers that don't ask for this (e.g. the interactive + // tool, which lets the user pan beyond the initial view) + uint8_t* cell_in_view = NULL; + if(restrict_mesh_azimuth) + { + cell_in_view = malloc((size_t)(2*render_radius_cells) * (size_t)(2*render_radius_cells)); + if(cell_in_view == NULL) + { + MSG("malloc(cell_in_view) failed"); + goto done; + } + + // Same formula as horizonator_move() uses to place the viewer + // within the loaded cell grid + const float viewer_cell_i = + (viewer_lon - ctx->dems.origin_dem_lon_lat[0]) * ctx->dems.cells_per_deg - + ctx->dems.origin_dem_cellij[0]; + const float viewer_cell_j = + (viewer_lat - ctx->dems.origin_dem_lon_lat[1]) * ctx->dems.cells_per_deg - + ctx->dems.origin_dem_cellij[1]; + const float cos_viewer_lat_local = cosf(viewer_lat * (float)M_PI/180.0f); + + const float az_center = (mesh_az_deg0 + mesh_az_deg1)/2.0f; + const float az_lo = mesh_az_deg0 - MESH_AZIMUTH_MARGIN_DEG; + const float az_hi = mesh_az_deg1 + MESH_AZIMUTH_MARGIN_DEG; + + const int W = 2*render_radius_cells; + for(int j=0; j= az_lo && az_deg <= az_hi); + } + + cell_in_view[j*W + i] = in_view ? 1 : 0; + } + } + // indices { GLuint indexBufID; @@ -502,10 +637,25 @@ bool horizonator_init( // output GLuint* indices = glMapBuffer(GL_ELEMENT_ARRAY_BUFFER, GL_WRITE_ONLY); int idx = 0; + const int W = 2*render_radius_cells; for( int j=0; j<(2*render_radius_cells-1); j++ ) { for( int i=0; i<(2*render_radius_cells-1); i++ ) { + if(cell_in_view != NULL) + { + // Skip this quad entirely unless at least one of its 4 + // corners is in view. This may keep a thin sliver of + // extra triangles right at the boundary; that's fine + bool any_in_view = + cell_in_view[(j+0)*W + (i+0)] || + cell_in_view[(j+1)*W + (i+1)] || + cell_in_view[(j+1)*W + (i+0)] || + cell_in_view[(j+0)*W + (i+1)]; + if(!any_in_view) + continue; + } + indices[idx++] = (j + 0)*(2*render_radius_cells) + (i + 0); indices[idx++] = (j + 1)*(2*render_radius_cells) + (i + 1); indices[idx++] = (j + 1)*(2*render_radius_cells) + (i + 0); @@ -517,7 +667,13 @@ bool horizonator_init( // output } int res = glUnmapBuffer(GL_ELEMENT_ARRAY_BUFFER); assert( res == GL_TRUE ); - assert(idx == ctx->Ntriangles*3); + assert(idx <= ctx->Ntriangles*3); + + // cell_in_view may have shrunk the actual triangle count below the + // worst-case estimate used to size the buffer above + ctx->Ntriangles = idx/3; + + free(cell_in_view); } // shaders @@ -614,6 +770,11 @@ bool horizonator_init( // output ctx->uniform_zfar = glGetUniformLocation(ctx->program, "zfar"); assert_opengl(); ctx->uniform_znear_color = glGetUniformLocation(ctx->program, "znear_color"); assert_opengl(); ctx->uniform_zfar_color = glGetUniformLocation(ctx->program, "zfar_color"); assert_opengl(); + ctx->uniform_curvature_scale = glGetUniformLocation(ctx->program, "curvature_scale"); assert_opengl(); + ctx->uniform_refraction_k = glGetUniformLocation(ctx->program, "refraction_k"); assert_opengl(); + ctx->uniform_shading_scale = glGetUniformLocation(ctx->program, "shading_scale"); assert_opengl(); + ctx->uniform_sun_dir = glGetUniformLocation(ctx->program, "sun_dir"); assert_opengl(); + ctx->uniform_materials_scale = glGetUniformLocation(ctx->program, "materials_scale"); assert_opengl(); #undef make_and_set_uniform // And I set the other uniforms @@ -621,6 +782,12 @@ bool horizonator_init( // output horizonator_set_zextents(ctx, HORIZONATOR_ZNEAR_DEFAULT, HORIZONATOR_ZFAR_DEFAULT, HORIZONATOR_ZNEAR_DEFAULT, HORIZONATOR_ZFAR_DEFAULT); + // Curvature correction is off by default: unchanged legacy behavior + horizonator_set_curvature(ctx, false, 0.13f); + // Slope shading is off by default: unchanged legacy behavior + horizonator_set_sun(ctx, false, 135.0f, 45.0f); + // Procedural materials are off by default: unchanged legacy behavior + horizonator_set_materials(ctx, false); } if(offscreen_width > 0) @@ -695,6 +862,14 @@ void horizonator_deinit( horizonator_context_t* ctx ) glutDestroyWindow(ctx->glut_window); ctx->glut_window = 0; } + + // Was leaked previously: the DEM files stay mmap-ed (and their pages + // resident, once touched by horizonator_dem_sample()) until this is + // called. This matters most for a process that calls + // horizonator_init()/horizonator_deinit() in a loop (e.g. to render + // several tiles one at a time): without this, memory use grows + // unboundedly across iterations + horizonator_dem_deinit(&ctx->dems); } bool horizonator_move(horizonator_context_t* ctx, @@ -893,6 +1068,65 @@ bool horizonator_set_zextents(horizonator_context_t* ctx, return true; } +bool horizonator_set_curvature(horizonator_context_t* ctx, + bool curvature_enabled, + float refraction_k) +{ + if(ctx->use_glut) + { + if(ctx->glut_window == 0) + return false; + glutSetWindow(ctx->glut_window); + } + + glUniform1f( ctx->uniform_curvature_scale, curvature_enabled ? 1.0f : 0.0f); assert_opengl(); + glUniform1f( ctx->uniform_refraction_k, refraction_k); assert_opengl(); + + return true; +} + +bool horizonator_set_sun(horizonator_context_t* ctx, + bool shading_enabled, + float sun_az_deg, + float sun_el_deg) +{ + if(ctx->use_glut) + { + if(ctx->glut_window == 0) + return false; + glutSetWindow(ctx->glut_window); + } + + // Same (east,north,height) convention as vertex.glsl: az=0 is North, + // az=90 is East (az_rad = atan(east,north) there) + const float az_rad = sun_az_deg * (float)M_PI/180.0f; + const float el_rad = sun_el_deg * (float)M_PI/180.0f; + const float cos_el = cos(el_rad); + const float east = cos_el * sin(az_rad); + const float north = cos_el * cos(az_rad); + const float height = sin(el_rad); + + glUniform1f( ctx->uniform_shading_scale, shading_enabled ? 1.0f : 0.0f); assert_opengl(); + glUniform3f( ctx->uniform_sun_dir, east, north, height); assert_opengl(); + + return true; +} + +bool horizonator_set_materials(horizonator_context_t* ctx, + bool materials_enabled) +{ + if(ctx->use_glut) + { + if(ctx->glut_window == 0) + return false; + glutSetWindow(ctx->glut_window); + } + + glUniform1f( ctx->uniform_materials_scale, materials_enabled ? 1.0f : 0.0f); assert_opengl(); + + return true; +} + bool horizonator_redraw(const horizonator_context_t* ctx) { if(ctx->use_glut) @@ -1119,7 +1353,10 @@ bool horizonator_project( // output double az_rad0, double az_rad1, int width, - int height) + int height, + + bool curvature_enabled, + double refraction_k) { const float Rearth = 6371000.0; @@ -1148,9 +1385,16 @@ bool horizonator_project( // output // The projection code is mostly lifted from vertex.glsl. Would be nice to // consolidate - const double h = ele - ele_viewer; const double distance_ne = sqrt(distance_sq_ne); - *range = sqrt(distance_sq_ne + h*h); + + // Same curvature-and-refraction correction as vertex.glsl. Must match, + // or this projection won't agree with the actual render + const double drop = + curvature_enabled ? + (1.0 - refraction_k) * distance_sq_ne / (2.0*Rearth) : 0.0; + const double h = (ele - ele_viewer) - drop; + + *range = sqrt(distance_sq_ne + h*h); const double aspect = (double)width / (double)height; const double el_ndc = atan2(h, distance_ne) * aspect * az_ndc_per_rad; diff --git a/horizonator-pywrap.c b/horizonator-pywrap.c index e50c934..ce661e2 100644 --- a/horizonator-pywrap.c +++ b/horizonator-pywrap.c @@ -109,6 +109,7 @@ py_horizonator_init(py_horizonator_t* self, PyObject* args, PyObject* kwargs) NULL, width, height, render_radius_cells, render_radius_m, + false, 0., 0., // no mesh azimuth restriction: render() may use any azimuth true, render_texture, SRTM1, dir_dems, dir_tiles, diff --git a/horizonator.cc b/horizonator.cc index eefd78a..8d26fd4 100644 --- a/horizonator.cc +++ b/horizonator.cc @@ -234,6 +234,7 @@ class GLWidget : public Fl_Gl_Window NULL, -1, -1, -1, zfar, + false, 0.f, 0.f, // no mesh azimuth restriction: I pan around freely false, render_texture, SRTM1, NULL,NULL, diff --git a/horizonator.h b/horizonator.h index db32c29..7230edf 100644 --- a/horizonator.h +++ b/horizonator.h @@ -33,6 +33,9 @@ typedef struct int32_t uniform_texturemap_dlat2; int32_t uniform_znear, uniform_zfar; int32_t uniform_znear_color, uniform_zfar_color; + int32_t uniform_curvature_scale, uniform_refraction_k; + int32_t uniform_shading_scale, uniform_sun_dir; + int32_t uniform_materials_scale; uint32_t program; @@ -81,6 +84,18 @@ static bool horizonator_context_isvalid(const horizonator_context_t* ctx) // SRTM1 selects between 1" SRTM and 3" SRTM. Currently every triangle is // rendered, so 1" SRTM tiles can easily overload the machine. Unless you need // the extra resolution, stick with 3" SRTM tiles for now +// +// By default the mesh covers the full circle of loaded DEM data (radius +// render_radius_cells/render_radius_m), even though a given render usually +// only looks at a fraction of that circle (see horizonator_pan_zoom()). If +// restrict_mesh_azimuth is true, the mesh is built only for the +// [mesh_az_deg0,mesh_az_deg1] azimuth wedge (plus a small margin), which can +// dramatically cut the triangle count for a narrow panorama. This trades +// away the ability to horizonator_pan_zoom() outside that wedge later +// without gaps in the mesh, so it should only be used when the caller knows +// it will never look outside that wedge (e.g. the "standalone" tool, which +// renders one fixed view). Leave false for callers that let the user pan +// around after loading (e.g. the interactive "horizonator" tool) bool horizonator_init( // output horizonator_context_t* ctx, @@ -95,6 +110,9 @@ bool horizonator_init( // output int render_radius_cells, // This should be given >0 float render_radius_m, // or this, but not both + bool restrict_mesh_azimuth, + float mesh_az_deg0, float mesh_az_deg1, // ignored unless restrict_mesh_azimuth + bool use_glut, bool render_texture, bool SRTM1, @@ -139,6 +157,37 @@ bool horizonator_set_zextents(horizonator_context_t* ctx, float znear, float zfar, float znear_color, float zfar_color); +// Enables/disables the Earth-curvature-and-refraction correction to the +// apparent elevation angle of rendered terrain. When curvature_enabled is +// false (the default set by horizonator_init()), rendering uses the +// original flat tangent-plane approximation, unchanged. refraction_k is +// the atmospheric refraction coefficient (~0.13 is a commonly-used value; +// see udeuschle.de); it is ignored when curvature_enabled is false +bool horizonator_set_curvature(horizonator_context_t* ctx, + bool curvature_enabled, + float refraction_k); + +// Enables/disables slope shading: terrain is darkened/lightened by a +// directional light, based on a smoothly-interpolated per-vertex surface +// normal (estimated from the DEM). When shading_enabled is false (the +// default set by horizonator_init()), rendering is unchanged (the plain +// distance-based grayscale). sun_az_deg (0=North, 90=East) and sun_el_deg +// (0=horizon, 90=straight up) give the direction TOWARDS the sun; both +// are ignored when shading_enabled is false +bool horizonator_set_sun(horizonator_context_t* ctx, + bool shading_enabled, + float sun_az_deg, + float sun_el_deg); + +// Enables/disables a first, purely procedural land-cover approximation +// (no aerial imagery or real land-cover data): each point is tinted by +// elevation (snow above a fixed snow line) and slope steepness (bare rock +// on steep terrain), forest below the tree line and alpine grass above it +// otherwise. When materials_enabled is false (the default set by +// horizonator_init()), rendering is unchanged +bool horizonator_set_materials(horizonator_context_t* ctx, + bool materials_enabled); + bool horizonator_redraw(const horizonator_context_t* ctx); // returns true if an intersection is found @@ -177,6 +226,11 @@ bool horizonator_x_from_az( // output double az_rad1, int width); +// curvature_enabled/refraction_k must match whatever was passed to +// horizonator_set_curvature() for the render being annotated -- otherwise +// the predicted screen position of (lat,lon,ele) will be off by the +// curvature drop (can be >1km of apparent height at 100+km), and this will +// incorrectly look occluded/unmatched bool horizonator_project( // output double* x, double* y, @@ -193,7 +247,10 @@ bool horizonator_project( // output double az_rad0, double az_rad1, int width, - int height); + int height, + + bool curvature_enabled, + double refraction_k); bool horizonator_unproject(// output float* lat, float* lon, diff --git a/query-peaks-from-osm.py b/query-peaks-from-osm.py index ae19653..6e5b152 100755 --- a/query-peaks-from-osm.py +++ b/query-peaks-from-osm.py @@ -42,10 +42,10 @@ def parse_args(): import requests import json -api = 'http://overpass-api.de/api/interpreter' +api = 'https://overpass-api.de/api/interpreter' query = f''' -[out:json]; +[out:json][timeout:60]; node ["natural" = "peak" ] @@ -55,8 +55,13 @@ def parse_args(): out; ''' +# The Overpass API rejects requests with no (or a generic) User-Agent, and +# plain http:// tends to time out; https:// with an identifying UA is +# reliable r = requests.post(api, - data = query) + data = query, + headers = {'User-Agent': 'horizonator query-peaks-from-osm.py'}, + timeout = 90) if r.status_code != 200: print(f"Overpass api error: {r.status_code} {r.reason=}", file = sys.stderr) diff --git a/standalone.c b/standalone.c index 487cda8..57515b6 100644 --- a/standalone.c +++ b/standalone.c @@ -12,8 +12,61 @@ #include "annotator.h" #include "util.h" +// Post-processing pass, working directly off the range (depth) buffer: +// darkens to a dark gray any pixel that sits at a depth discontinuity -- +// either terrain against the sky, or one terrain surface abruptly giving +// way to another, much nearer or farther one behind a gap (a col, a +// notch...). This traces a crisp line along the crest of each visible +// ridge/mountain range "layer", the way classic drawn panoramas (e.g. +// udeuschle.de) do it. Only checks vertically (row to row): ridge crests +// are near-horizontal features in a panorama, so this is where the sharp +// depth jumps are. +// +// The jump is judged as an absolute distance, in meters: at a fixed +// threshold, grazing, near-horizontal sightlines (typically towards the +// bottom of the frame, at low elevation angles) produce many more of +// these lines than steeper sightlines towards distant peaks, since a +// one-pixel vertical step there corresponds to a much bigger jump in true +// distance. This is no different than in a real drawn panorama: nearby +// terrain, seen at a grazing angle, really does present many more +// distinct ridgelets per degree of elevation than distant terrain does +static void draw_ridge_outlines(uint8_t* image, const float* ranges, + int width, int height, + float threshold_m, + float ridge_line_gray) +{ + const uint8_t gray_u8 = (uint8_t)(ridge_line_gray*255.0f + 0.5f); + + for(int x=0; x threshold_m); // one surface replaced by another + + if(edge) + { + uint8_t* px = &image[((size_t)y*width + x)*3]; + px[0] = px[1] = px[2] = gray_u8; + } + } +} + static bool glut_loop( bool render_texture, bool SRTM1, float viewer_lat, float viewer_lon, + float viewer_height_m, + bool curvature_enabled, float refraction_k, + bool shading_enabled, float sun_az_deg, float sun_el_deg, + bool materials_enabled, + bool restrict_mesh_azimuth, // Bounds of the view. We expect az_deg1 > az_deg0. The azimuth // edges lie at the edges of the image. So for an image that's @@ -35,11 +88,13 @@ static bool glut_loop( bool render_texture, bool SRTM1, { horizonator_context_t ctx; + float viewer_z = -1.0f; if( !horizonator_init( &ctx, viewer_lat, viewer_lon, - NULL, + &viewer_z, -1, -1, -1, zfar, + restrict_mesh_azimuth, az_deg0, az_deg1, true, render_texture, SRTM1, dir_dems, @@ -49,6 +104,25 @@ static bool glut_loop( bool render_texture, bool SRTM1, allow_downloads) ) return false; + if(viewer_height_m != 0.0f) + { + // horizonator_init() auto-selected a viewer_z sitting on the DEM + // ground surface. I add the height of the observer above that + // ground (e.g. the floor of an apartment building) and re-apply it. + viewer_z += viewer_height_m; + if(!horizonator_move(&ctx, &viewer_z, viewer_lat, viewer_lon)) + return false; + } + + if(!horizonator_set_curvature(&ctx, curvature_enabled, refraction_k)) + return false; + + if(!horizonator_set_sun(&ctx, shading_enabled, sun_az_deg, sun_el_deg)) + return false; + + if(!horizonator_set_materials(&ctx, materials_enabled)) + return false; + if(!horizonator_set_zextents(&ctx, znear, zfar, znear_color, zfar_color)) return false; @@ -114,59 +188,202 @@ int main(int argc, char* argv[]) { const char* usage = "%s [--width WIDTH_PIXELS] [--height HEIGHT_PIXELS]\n" - " [--image OUT.png|OUT.pdf|OUT.svg]\n" + " [--image OUT.png|OUT.pdf|OUT.svg] [--cut-off-bottom-px N]\n" " [--texture] [--SRTM1]\n" " [--allow-tile-downloads]\n" - " [--znear ZNEAR]\n" - " [--zfar ZFAR]\n" - " [--znear-color ZNEARCOLOR]\n" - " [--zfar-color ZFARCOLOR]\n" - " [--dirdems DIRECTORY]\n" - " [--dirtiles DIRECTORY]\n" - " [--tiles NAME=FMT]\n" + " [--viewer-height METERS]\n" + " [--curvature] [--refraction-k K]\n" + " [--shading] [--sun-azimuth DEG] [--sun-elevation DEG]\n" + " [--materials]\n" + " [--no-restrict-mesh-azimuth]\n" + " [--no-ridge-lines] [--ridge-line-threshold METERS] [--ridge-line-gray FRACTION]\n" + " [--label-font-size-pt POINTS]\n" + " [--list-visible-peaks OUT.txt]\n" + " [--znear ZNEAR] [--zfar ZFAR] [--znear-color ZNEARCOLOR] [--zfar-color ZFARCOLOR]\n" + " [--dirdems DIRECTORY] [--dirtiles DIRECTORY] [--tiles NAME=FMT]\n" " LAT LON AZ_CENTER_DEG AZ_RADIUS_DEG\n" "\n" + "=== Basic operation ===\n" + "\n" "By default, we render to a window. If --width is given, we render\n" - "to an image (--image) instead.\n" - "--height applies only if --width is given, and is optional; a reasonable\n" - "field-of-view will be assumed if --height is omitted." - "If we're annotating a --image, you can pass --cut-off-bottom-px to cut off\n" - "the given number of pixels from the bottom. This is a workaround for the\n" - "uneven edges of the render at the bottom\n" + "to an image (--image) instead. --height applies only if --width is\n" + "given, and is optional; a reasonable field-of-view is assumed if\n" + "--height is omitted (see FIELD OF VIEW below).\n" + "\n" + "LAT, LON are the coordinates of the viewer, in degrees. AZ_CENTER_DEG\n" + "and AZ_RADIUS_DEG give the azimuth window we look at: from\n" + "AZ_CENTER_DEG-AZ_RADIUS_DEG to AZ_CENTER_DEG+AZ_RADIUS_DEG (0 = North,\n" + "90 = East). When plotting to a window, these describe the azimuth\n" + "bounds of the VIEWPORT; when rendering to an image, the centers of\n" + "the first and last pixels (slightly narrower than the viewport: one\n" + "extra half-pixel on each side).\n" + "\n" + "The image filename MUST be a .png file (the render alone is written)\n" + "OR a .pdf or .svg file (the render, annotated with named peaks and\n" + "bearing markers, is written -- see PEAK LABELS below). If annotating,\n" + "--cut-off-bottom-px N discards the bottom N pixels of the render; a\n" + "workaround for the uneven edges some renders have at the bottom.\n" + "\n" + "=== Field of view: what to render, and how much of it fills the frame ===\n" "\n" - "The image filename MUST be a .png file (the render will be written)\n" - "OR a .pdf or .svg file (the annotated render will be written)\n" + "--znear/--zfar (meters) bound how much of the terrain is rendered:\n" + "closer than --znear or farther than --zfar is not drawn at all. Both\n" + "have reasonable defaults and may be omitted. Raising --znear is a\n" + "good way to skip an uninteresting nearby plain, dedicating the whole\n" + "frame to the distant relief that's actually worth looking at.\n" "\n" - "When plotting to a window, AZ_..._DEG refers to the azimuth bounds of the\n" - "VIEWPORT. When rendering to an image, to the\n" - "centers of the first and last pixels. This is slightly smaller\n" - "than the whole viewport: there's one extra pixel on each side\n" + "The image's vertical field of view is NOT set directly: it falls out\n" + "of --width, --height and AZ_RADIUS_DEG, because the same angular\n" + "resolution (degrees/pixel) is used both horizontally and vertically:\n" + " vertical_FOV_deg = 2*AZ_RADIUS_DEG * height/width\n" + "A real mountain panorama usually only occupies a few degrees of\n" + "elevation above/below the horizontal, so with the default --height\n" + "(a 20deg half-angle assumed when omitted) most of the frame ends up\n" + "as empty sky. Passing an explicit, small --height (to get a vertical\n" + "FOV of, say, 5-10deg) both crops that away AND increases the\n" + "effective resolution devoted to the relief that matters -- this is\n" + "the single biggest lever for a sharp-looking render. The same logic\n" + "applies to AZ_RADIUS_DEG: narrowing it, at a fixed --width, increases\n" + "the angular resolution (degrees/pixel) of that slice.\n" "\n" - "The extents of ranges that we render are given by --znear and --zfar,\n" - "in meters. Anything closer than --znear and further out than --zfar will\n" - "not be rendered. The color-coding extents are given by --znear-color and\n" - "--zfar-color. Anything at or closer than --znear-color will be rendered\n" - "as black, and anything at or further than --zfar-color will be rendered\n" - "as red. Of --zfar-color/--znear-color are omitted, they take the same value\n" - "as --zfar/--znear. --znear and --zfar have reasonable defaults, and may also\n" - "be omitted\n" + "=== Earth curvature and atmospheric refraction ===\n" "\n" - "By default we colorcode the renders by range. If --texture, we\n" - "use a set of image tiles to texture the render instead\n" + "By default we render assuming a flat tangent plane, ignoring the\n" + "curvature of the Earth -- a reasonable approximation at short range,\n" + "but at 100+km (distant high peaks) the apparent elevation angle can\n" + "be off by upwards of a kilometer of apparent height. Pass --curvature\n" + "to correct for it:\n" + " drop = (1-k) * distance^2 / (2*Rearth)\n" + "subtracted from the apparent height of each rendered point.\n" + "--refraction-k sets the refraction coefficient k (default 0.13, a\n" + "commonly-used value; see udeuschle.de). It is ignored unless\n" + "--curvature is also given. k=0 is pure geometric curvature, with no\n" + "refraction compensation (peaks look their lowest); increasing k\n" + "raises them back up somewhat, closer to the flat-plane render.\n" "\n" - "By default we use 3\" SRTM data. Currently every triangle in the grid is\n" - "rendered. This is inefficient, but the higher-resolution 1\" SRTM tiles\n" - "would make it use 9 times more memory and computational resources, so\n" - "sticking with the lower-resolution 3\" SRTM data is recommended for now.\n" + "=== Slope shading ===\n" + "\n" + "By default terrain is colored purely by distance (see COLOR below),\n" + "with no directional lighting. Pass --shading to darken/lighten the\n" + "relief by a directional light, based on a smoothly-interpolated\n" + "per-vertex surface normal (estimated from the DEM), giving the\n" + "relief a 3D appearance. --sun-azimuth (default 135, i.e. SE) and\n" + "--sun-elevation (default 45deg above the horizon) set the direction\n" + "the light comes from; both are ignored unless --shading is given.\n" + "A slope facing away from the sun is dimmed, never made fully black.\n" + "\n" + "=== Materials ===\n" + "\n" + "By default terrain isn't tinted by land cover, only by distance and\n" + "(if --shading) slope lighting. Pass --materials for a first, purely\n" + "procedural approximation: each point is classified by elevation\n" + "(snow above a fixed snow line) and slope steepness (bare rock on\n" + "steep terrain), with forest below the tree line and alpine grass\n" + "above it otherwise. This uses no aerial imagery or real land-cover\n" + "data, just the same DEM already being rendered -- a quick\n" + "placeholder, not a substitute for real land-cover texturing.\n" + "\n" + "=== Color and ridge outlines ===\n" + "\n" + "By default the render is color-coded by range (atmospheric\n" + "perspective): near terrain is drawn in a dark neutral gray, far\n" + "terrain fades towards a pale blue-gray, closer to the white\n" + "background, the way haze tints distant relief in a real photo.\n" + "--znear-color/--zfar-color set the distance extents used for this\n" + "(in meters); if omitted, they take the same values as --znear/--zfar.\n" + "\n" + "On top of that, the crest of each visible ridge/mountain-range\n" + "\"layer\" is outlined, the way classic drawn panoramas (e.g.\n" + "udeuschle.de) do it: a post-process on the range image darkens any\n" + "pixel that sits at a sharp depth discontinuity -- terrain against\n" + "the sky, or one surface abruptly giving way to a much nearer/farther\n" + "one behind a gap. Pass --no-ridge-lines to disable this.\n" + "--ridge-line-threshold (meters, default 500) sets how big a range\n" + "jump between vertically-adjacent pixels counts as a new layer: lower\n" + "values pick out more (and finer) ridgelets, especially at low\n" + "elevation angles / grazing sightlines near the bottom of the frame,\n" + "where a one-pixel step naturally corresponds to a much bigger jump\n" + "in true distance than the same step towards a distant peak; this is\n" + "expected, not a bug, and mirrors what a real drawn panorama looks\n" + "like up close. --ridge-line-gray (0-1, default 0.15) sets how dark\n" + "the outline is: 0 is black, 1 is white. Keep it darker than the\n" + "nearest terrain gray so it stays visible against it.\n" + "\n" + "If --texture, we use a set of image tiles to texture the render\n" + "instead of any of the above color-coding.\n" + "\n" + "=== Peak labels ===\n" + "\n" + "See query-peaks-from-osm.py to generate the compiled-in list of\n" + "named peaks (a poi_t[] #included by this file -- edit the #include\n" + "line to point at your own generated file). A peak counts as visible\n" + "if it's (a) actually found in the rendered range image (not\n" + "occluded by closer terrain), (b) farther than 500m, and (c) no\n" + "farther than --zfar. --viewer-height and --curvature/--refraction-k\n" + "affect where a peak is predicted to land on screen; this tool\n" + "always uses the same values for the render and for the labelling,\n" + "so there's nothing extra to keep in sync here.\n" + "\n" + "--list-visible-peaks OUT.txt writes the visible peaks (independent\n" + "of --image; works even without one) to a tab-separated text file:\n" + "name, lat, lon, ele_m, range_m, one per line, in no particular\n" + "order. Peaks with no OSM name fall back to their altitude as a\n" + "name (e.g. \"1583.0\") -- filter those out downstream by checking\n" + "for at least one letter, if wanted.\n" + "\n" + "With a .pdf/.svg --image, visible peaks are labelled on the\n" + "render: a black name in vertical text (read by tilting your head\n" + "to the left: first letter at the bottom, last letter at the top),\n" + "connected to its peak by a thin almond-green leader line. Each\n" + "label's last (top) character sits 30px below the top of the\n" + "canvas, regardless of how far its peak actually is -- so the\n" + "leader line length alone shows how far up the label had to reach\n" + "in the sky. Before drawing, peaks whose screen-x positions are\n" + "closer than --label-font-size-pt*1.5 are considered to conflict;\n" + "within each conflicting group, only the highest-elevation peak\n" + "gets a label, the others are dropped entirely (no line, no name).\n" + "--label-font-size-pt sets the text height in real typographic\n" + "points, i.e. as it will appear in the output PDF/SVG page (default\n" + "12); raise it for a legible render at a large --width, or when the\n" + "page will be viewed/printed at less than 100% zoom -- this also\n" + "widens the conflict-detection gap, so fewer, larger labels survive.\n" + "\n" + "=== Viewer position ===\n" + "\n" + "--viewer-height adds the given number of meters to the viewer\n" + "elevation sampled from the DEM at LAT,LON (e.g. the height of an\n" + "apartment floor above street level). Defaults to 0.\n" + "\n" + "=== Performance ===\n" + "\n" + "By default we use 3\" SRTM data, and only mesh the AZ_CENTER_DEG +-\n" + "AZ_RADIUS_DEG wedge that's actually being rendered (plus a small\n" + "margin), rather than the full circle of loaded DEM data (radius\n" + "--zfar) -- this cuts the triangle count, and so the render time,\n" + "without changing the output, since this tool always renders the\n" + "same wedge it was given on the commandline. Pass\n" + "--no-restrict-mesh-azimuth to mesh the full circle instead (slower;\n" + "only useful for debugging).\n" + "\n" + "The higher-resolution 1\" SRTM tiles can be selected with --SRTM1;\n" + "this gives a 9x increase in triangle count (currently every\n" + "triangle in the mesh is rendered, however far away and however tiny\n" + "on screen), so it can easily overload the machine. Stick with the\n" + "default 3\" data unless you need the extra resolution for very\n" + "close terrain.\n" + "\n" + "=== Data sources ===\n" "\n" "The DEMs are in the directory given by --dirdems, or in\n" "~/.horizonator/DEMs_SRTM3/ (or DEMs_SRTM1) if omitted.\n" "\n" "The tiles are in the directory given by --dirtiles, or in\n" - "~/.horizonator/tiles if omitted. This is the BASE directory for ALL the\n" - "available tile sets. By default we use the OSM mapnik tiles. To specify\n" - "different tiles, pass '--tiles NAME=FMT'. Where NAME is the identifier of\n" - "this set and FMT is the URL format string to use for this set\n"; + "~/.horizonator/tiles if omitted. This is the BASE directory for ALL\n" + "the available tile sets. By default we use the OSM mapnik tiles. To\n" + "specify different tiles, pass '--tiles NAME=FMT', where NAME is the\n" + "identifier of this set and FMT is the URL format string to use for\n" + "this set. --allow-tile-downloads lets missing tiles be fetched over\n" + "the network; without it, a missing tile is just left blank.\n"; struct option opts[] = { { "width", required_argument, NULL, 'w' }, @@ -179,6 +396,19 @@ int main(int argc, char* argv[]) { "texture", no_argument, NULL, 'T' }, { "SRTM1", no_argument, NULL, 'S' }, { "allow-tile-downloads",no_argument, NULL, 'a' }, + { "viewer-height", required_argument, NULL, 'V' }, + { "curvature", no_argument, NULL, 'C' }, + { "refraction-k", required_argument, NULL, 'k' }, + { "shading", no_argument, NULL, 's' }, + { "sun-azimuth", required_argument, NULL, 'A' }, + { "sun-elevation", required_argument, NULL, 'E' }, + { "materials", no_argument, NULL, 'm' }, + { "no-restrict-mesh-azimuth", no_argument, NULL, 'M' }, + { "no-ridge-lines", no_argument, NULL, 'R' }, + { "ridge-line-threshold", required_argument, NULL, 'r' }, + { "ridge-line-gray", required_argument, NULL, 'g' }, + { "label-font-size-pt", required_argument, NULL, 'F' }, + { "list-visible-peaks", required_argument, NULL, 'L' }, { "znear", required_argument, NULL, '1' }, { "zfar", required_argument, NULL, '2' }, { "znear-color", required_argument, NULL, '3' }, @@ -191,6 +421,7 @@ int main(int argc, char* argv[]) int height = 0; int cut_off_bottom_px = 0; const char* filename_image = NULL; + const char* filename_visible_peaks = NULL; const char* dir_dems = NULL; const char* dir_tiles = NULL; const char* tiles_name = NULL; @@ -198,6 +429,18 @@ int main(int argc, char* argv[]) bool render_texture = false; bool SRTM1 = false; bool allow_downloads = false; + float viewer_height_m = 0.0f; + bool curvature_enabled = false; + float refraction_k = 0.13f; + bool shading_enabled = false; + float sun_az_deg = 135.0f; + float sun_el_deg = 45.0f; + bool materials_enabled = false; + bool restrict_mesh_azimuth = true; + bool ridge_lines = true; + float ridge_line_threshold_m = 500.0f; + float ridge_line_gray = 0.15f; + float label_font_size_pt = 12.0f; float znear = HORIZONATOR_ZNEAR_DEFAULT; float zfar = HORIZONATOR_ZFAR_DEFAULT; @@ -315,6 +558,58 @@ int main(int argc, char* argv[]) allow_downloads = true; break; + case 'V': + viewer_height_m = (float)atof(optarg); + break; + + case 'C': + curvature_enabled = true; + break; + + case 'k': + refraction_k = (float)atof(optarg); + break; + + case 's': + shading_enabled = true; + break; + + case 'A': + sun_az_deg = (float)atof(optarg); + break; + + case 'E': + sun_el_deg = (float)atof(optarg); + break; + + case 'm': + materials_enabled = true; + break; + + case 'M': + restrict_mesh_azimuth = false; + break; + + case 'R': + ridge_lines = false; + break; + + case 'r': + ridge_line_threshold_m = (float)atof(optarg); + break; + + case 'g': + ridge_line_gray = (float)atof(optarg); + break; + + case 'F': + label_font_size_pt = (float)atof(optarg); + break; + + case 'L': + filename_visible_peaks = optarg; + break; + case '?': fprintf(stderr, "Unknown option\n\n"); fprintf(stderr, usage, argv[0]); @@ -333,15 +628,15 @@ int main(int argc, char* argv[]) if(znear_color < 0.f) znear_color = znear; if(zfar_color < 0.f) zfar_color = zfar; - if(width > 0 && filename_image == NULL) + if(width > 0 && filename_image == NULL && filename_visible_peaks == NULL) { - fprintf(stderr, "--width makes sense only with --image\n\n"); + fprintf(stderr, "--width makes sense only with --image or --list-visible-peaks\n\n"); fprintf(stderr, usage, argv[0]); return 1; } - if(width <= 0 && filename_image != NULL) + if(width <= 0 && (filename_image != NULL || filename_visible_peaks != NULL)) { - fprintf(stderr, "--width required if --image\n\n"); + fprintf(stderr, "--width required if --image or --list-visible-peaks\n\n"); fprintf(stderr, usage, argv[0]); return 1; } @@ -370,10 +665,15 @@ int main(int argc, char* argv[]) } - if(filename_image == NULL) + if(filename_image == NULL && filename_visible_peaks == NULL) { glut_loop(render_texture, SRTM1, lat, lon, + viewer_height_m, + curvature_enabled, refraction_k, + shading_enabled, sun_az_deg, sun_el_deg, + materials_enabled, + restrict_mesh_azimuth, az_center_deg-az_radius_deg, az_center_deg+az_radius_deg, znear,zfar,znear_color,zfar_color, @@ -383,9 +683,13 @@ int main(int argc, char* argv[]) return 0; } - const int strlen_filename_image = strlen(filename_image); + // filename_image can legitimately be NULL here now (--list-visible-peaks + // without --image); strlen(NULL) is undefined behavior, so this has to + // stay inside the NULL check, unlike before + int strlen_filename_image = 0; if(filename_image != NULL) { + strlen_filename_image = strlen(filename_image); if(!(strlen_filename_image >= 5 && (0 == strcasecmp(".png", &filename_image[strlen_filename_image-4]) || 0 == strcasecmp(".pdf", &filename_image[strlen_filename_image-4]) || @@ -413,7 +717,7 @@ int main(int argc, char* argv[]) uint8_t* pool = NULL; char* image; float* ranges; - if(filename_image != NULL) + if(filename_image != NULL || filename_visible_peaks != NULL) { // rgb for the image and float for the depth pool = malloc( width*height * (3 + sizeof(float)) ); @@ -435,6 +739,9 @@ int main(int argc, char* argv[]) &viewer_z, width, height, -1, zfar, + restrict_mesh_azimuth, + az_center_deg-az_radius_deg, + az_center_deg+az_radius_deg, true, render_texture, SRTM1, dir_dems, dir_tiles, @@ -445,6 +752,28 @@ int main(int argc, char* argv[]) return false; } + if(viewer_height_m != 0.0f) + { + // viewer_z was auto-selected by horizonator_init() to sit on the DEM + // ground surface. I add the height of the observer above that + // ground (e.g. the floor of an apartment building) and re-apply it. + viewer_z += viewer_height_m; + if(!horizonator_move(&ctx, &viewer_z, lat, lon)) + { + fprintf(stderr, "horizonator_move() failed\n"); + return false; + } + } + + if(!horizonator_set_curvature(&ctx, curvature_enabled, refraction_k)) + return false; + + if(!horizonator_set_sun(&ctx, shading_enabled, sun_az_deg, sun_el_deg)) + return false; + + if(!horizonator_set_materials(&ctx, materials_enabled)) + return false; + if(!horizonator_set_zextents(&ctx, znear, zfar, znear_color, zfar_color)) return false; @@ -463,6 +792,48 @@ int main(int argc, char* argv[]) return 1; } + if(ridge_lines) + draw_ridge_outlines((uint8_t*)image, ranges, width, height, + ridge_line_threshold_m, ridge_line_gray); + + // Shared by --image (.pdf/.svg) and --list-visible-peaks + poi_t pois[] = { +// ./query-peaks-from-osm.py 45.77294 4.82993 200000 > lyon-peaks.h +#include "socal-peaks.h" + }; + const int N_pois = (int)(sizeof(pois) / sizeof(pois[0])); + + if(filename_visible_peaks != NULL) + { + visible_poi_t visible[N_pois]; + int Nvisible = find_visible_pois(visible, + ranges, width, height, cut_off_bottom_px, + pois, N_pois, + lat, lon, + az_center_deg-az_radius_deg, + az_center_deg+az_radius_deg, + viewer_z, + curvature_enabled, refraction_k, + zfar); + + FILE* fp = fopen(filename_visible_peaks, "w"); + if(fp == NULL) + { + fprintf(stderr, "Couldn't open '%s' for writing\n", filename_visible_peaks); + return 1; + } + // namelatlonele_mrange_m, one visible peak per + // line, unsorted. See cluster-visible-peaks.py for a companion + // script that reads this format + for(int i=0; iname, poi->lat, poi->lon, poi->ele_m, visible[i].range); + } + fclose(fp); + } + if(filename_image != NULL) { if(0 == strcasecmp(".png", &filename_image[strlen_filename_image-4])) @@ -490,23 +861,20 @@ int main(int argc, char* argv[]) else { // pdf file is requested. I write an annotated pdf - poi_t pois[] = { -// ./query-peaks-from-osm.py 34. -118 100000 > socal-peaks.h -#include "socal-peaks.h" - }; - const int N_pois = (int)(sizeof(pois) / sizeof(pois[0])); - annotate(filename_image, (uint8_t*)image, ranges, width, height, cut_off_bottom_px, pois, N_pois, lat, lon, az_center_deg-az_radius_deg, az_center_deg+az_radius_deg, - viewer_z); + viewer_z, + curvature_enabled, refraction_k, + zfar, + label_font_size_pt); } - - free(pool); } + free(pool); + return 0; } diff --git a/vertex.glsl b/vertex.glsl index 8604204..47c2cac 100644 --- a/vertex.glsl +++ b/vertex.glsl @@ -3,6 +3,11 @@ #version 420 layout (location = 0) in vec3 vertex; +// Per-vertex normal (world-space east/north/height frame), estimated on +// the CPU by finite differences over the DEM (see horizonator-lib.c). Used +// for smooth slope shading: interpolated by the rasterizer across each +// triangle, so shading is continuous across triangle edges +layout (location = 1) in vec3 normal_attr; // We receive these from the CPU code uniform float viewer_cell_i, viewer_cell_j; @@ -23,10 +28,31 @@ uniform int osmtile_lowestX, osmtile_lowestY; uniform float znear, zfar; uniform float znear_color, zfar_color; +// Earth-curvature-and-refraction correction. curvature_scale is 0.0 +// (disabled: legacy flat tangent-plane rendering) or 1.0 (enabled). +// refraction_k is the atmospheric refraction coefficient (ignored if +// curvature_scale == 0) +uniform float curvature_scale; +uniform float refraction_k; + // We send these to the fragment shader -out vec3 rgb; +// +// atmo_t is 0 at znear_color and 1 at zfar_color: how far towards the +// atmospheric haze color this point should be blended, in fragment.glsl +// (which also knows about material colors, and needs this to decide the +// near/far blend from a common starting point -- see COLOR_NEAR_DEFAULT +// there) +out float atmo_t; out vec2 tex; +// Passed through to the geometry/fragment shaders for smooth slope shading +out vec3 normal; + +// Raw DEM elevation (meters above sea level, NOT relative to the viewer), +// for the procedural material classification (snow line etc.) in +// fragment.glsl +out float elevation_m; + const float Rearth = 6371000.0; const float pi = 3.14159265358979; @@ -128,9 +154,25 @@ void main(void) vec2 en = vec2( (i - viewer_cell_i) * DEG_PER_CELL * Rearth * pi/180. * cos_viewer_lat, (j - viewer_cell_j) * DEG_PER_CELL * Rearth * pi/180. ); - vec3 enh = vec3( en.x, en.y, vertex.z - viewer_z ); distance_ne = length(en); + + // A target at horizontal distance distance_ne appears lower than + // this flat-plane geometry predicts, because it sits behind the + // curve of the Earth; atmospheric refraction partially + // compensates by bending the light ray back down. The standard + // approximation for this net apparent drop is + // drop = (1-k) * d^2 / (2*Rearth) + // (k=0.13 is a commonly-used refraction coefficient; see + // udeuschle.de). curvature_scale==0 makes this vanish, giving + // back the original flat-plane behavior + float drop = curvature_scale * (1.0 - refraction_k) * + distance_ne*distance_ne / (2.0*Rearth); + + vec3 enh = vec3( en.x, en.y, vertex.z - viewer_z - drop ); + normal = normal_attr; + elevation_m = vertex.z; + float az_rad = atan(en.x, en.y); // az = 0: North @@ -156,8 +198,6 @@ void main(void) 1.0 ); } - rgb.r = max(min((distance_ne - znear_color) / (zfar_color - znear_color), - 1.0), 0.0); - rgb.g = 0.; - rgb.b = 0.; + atmo_t = clamp((distance_ne - znear_color) / (zfar_color - znear_color), + 0.0, 1.0); }