/*-------------------------------------------------------------------- * $Id: grdgradient.c 12822 2014-01-31 23:39:56Z remko $ * * Copyright (c) 1991-2014 by P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis and F. Wobbe * See LICENSE.TXT file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation; version 3 or any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * Brief synopsis: grdgradient read a grid file and compute gradient in azim direction: * azim = azimuth clockwise from north in degrees. * gradient = -[(dz/dx)sin(azim) + (dz/dy)cos(azim)]. * the expression in [] is the correct gradient. We take * -[] in order that data which goes DOWNHILL in the * azim direction will give a positive value; this is * for image shading purposes. Note for -D we instead use * the common convention of reporting the UPHILL direction. * * Author: W.H.F. Smith * Date: 1-JAN-2010 * Version: 5 API */ #define THIS_MODULE_NAME "grdgradient" #define THIS_MODULE_LIB "core" #define THIS_MODULE_PURPOSE "Compute directional gradients from a grid" #include "gmt_dev.h" #define GMT_PROG_OPTIONS "-RVfn" struct GRDGRADIENT_CTRL { struct In { bool active; char *file; } In; struct A { /* -A[/] */ bool active; bool two; double azimuth[2]; } A; struct D { /* -D[a][o][n] */ bool active; unsigned int mode; } D; struct E { /* -E[s|p]/ */ bool active; unsigned int mode; double azimuth, elevation; double ambient, diffuse, specular, shine; } E; struct G { /* -G */ bool active; char *file; } G; struct N { /* -N[t_or_e][[/[/]]] */ bool active; unsigned int mode; /* 1 = atan, 2 = exp */ double norm, sigma, offset; } N; struct S { /* -S */ bool active; char *file; } S; }; void *New_grdgradient_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ struct GRDGRADIENT_CTRL *C; C = GMT_memory (GMT, NULL, 1, struct GRDGRADIENT_CTRL); /* Initialize values whose defaults are not 0/false/NULL */ C->E.ambient = 0.55; C->E.diffuse = 0.6; C->E.specular = 0.4; C->E.shine = 10; C->N.norm = 1.0; return (C); } void Free_grdgradient_Ctrl (struct GMT_CTRL *GMT, struct GRDGRADIENT_CTRL *C) { /* Deallocate control structure */ if (!C) return; if (C->In.file) free (C->In.file); if (C->G.file) free (C->G.file); if (C->S.file) free (C->S.file); GMT_free (GMT, C); } double specular (double nx, double ny, double nz, double *s) { /* SPECULAR Specular reflectance. R = SPECULAR(Nx,Ny,Nz,S,V) returns the reflectance of a surface with normal vector components [Nx,Ny,Nz]. S and V specify the direction to the light source and to the viewer, respectively. For the time beeing I'm using V = [azim elev] = [0 90] so the following V[0] = sind(V[0])*cosd(V[1]); V[1] = -cosd(V[0])*cosd(V[1]); V[2] = sind(V[1]); Reduces to V[0] = 0; V[1] = 0; V[2] = 1 */ /*r = MAX(0,2*(s[0]*nx+s[1]*ny+s[2]*nz).*(v[0]*nx+v[1]*ny+v[2]*nz) - (v'*s)*ones(m,n)); */ return (MAX(0, 2 * (s[0]*nx + s[1]*ny + s[2]*nz) * nz - s[2])); } int GMT_grdgradient_usage (struct GMTAPI_CTRL *API, int level) { GMT_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); GMT_Message (API, GMT_TIME_NONE, "usage: grdgradient -G [-A[/]] [-D[a][o][n]]\n"); GMT_Message (API, GMT_TIME_NONE, "\t[-E[s|p|m]/[////]]\n"); GMT_Message (API, GMT_TIME_NONE, "\t[-N[t|e][[/[/]]]] [%s]\n\t[-S] [%s] [-fg] [%s]\n\n", GMT_Rgeo_OPT, GMT_V_OPT, GMT_n_OPT); if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Message (API, GMT_TIME_NONE, "\t is name of input grid file.\n"); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Message (API, GMT_TIME_NONE, "\t-A Set azimuth (0-360 CW from North (+y)) for directional derivatives.\n"); GMT_Message (API, GMT_TIME_NONE, "\t -A/ will compute two directions and save the one larger in magnitude.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-D Find the direction of the vector grad z (up-slope direction).\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append c to get cartesian angle (0-360 CCW from East (+x)) [Default: azimuth].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append o to get bidirectional orientations [0-180] rather than directions [0-360].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append n to add 90 degrees to the values from c or o.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-E Compute Lambertian radiance appropriate to use with grdimage/grdview.\n"); GMT_Message (API, GMT_TIME_NONE, "\t -E sets azimuth and elevation of light vector.\n"); GMT_Message (API, GMT_TIME_NONE, "\t -E///// sets azim, elev and\n"); GMT_Message (API, GMT_TIME_NONE, "\t other parameters that control the reflectance properties of the surface\n"); GMT_Message (API, GMT_TIME_NONE, "\t [Default values are: 0.55/0.6/0.4/10].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Specify '=' to get the default value (e.g., -E60/30/=/0.5).\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append s to use a simpler Lambertian algorithm (note that with this form\n"); GMT_Message (API, GMT_TIME_NONE, "\t you only have to provide the azimuth and elevation parameters).\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append p to use the Peucker piecewise linear approximation (simpler but faster algorithm).\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append m to use another algorithm that gives results close to ESRI's 'hillshade' but faster\n"); GMT_Message (API, GMT_TIME_NONE, "\t Note that in this case the azimuth and elevation are hardwired to 315 and 45 degrees.\n"); GMT_Message (API, GMT_TIME_NONE, "\t This means that even if you provide other values they will be ignored.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-G Output file for results from -A or -D.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-N Normalize gradients so that max |grad| = [1.0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Nt will make atan transform, then scale to [1.0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Ne will make exp transform, then scale to [1.0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t -Nt/[/] or -Ne/[/] sets sigma\n"); GMT_Message (API, GMT_TIME_NONE, "\t (and offset) for transform. [sigma, offset estimated from data].\n"); GMT_Option (API, "R"); GMT_Message (API, GMT_TIME_NONE, "\t-S Output file for |grad z|; requires -D.\n"); GMT_Option (API, "V"); GMT_Message (API, GMT_TIME_NONE, "\t-fg Convert geographic grids to meters using a \"Flat Earth\" approximation.\n"); GMT_Option (API, "n,."); return (EXIT_FAILURE); } int GMT_grdgradient_parse (struct GMT_CTRL *GMT, struct GRDGRADIENT_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to grdgradient and sets parameters in Ctrl. * Note Ctrl has already been initialized and non-zero default values set. * Any GMT common options will override values set previously by other commands. * It also replaces any file names specified as input or output with the data ID * returned when registering these sources/destinations with the API. */ unsigned int n_errors = 0, n_files = 0, j, entry, pos; int n_opt_args = 0; char ptr[GMT_BUFSIZ]; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; for (opt = options; opt; opt = opt->next) { /* Process all the options given */ switch (opt->option) { case '<': /* Input file (only one is accepted) */ if (n_files++ > 0) break; if ((Ctrl->In.active = GMT_check_filearg (GMT, '<', opt->arg, GMT_IN))) Ctrl->In.file = strdup (opt->arg); else n_errors++; break; /* Processes program-specific parameters */ case 'A': /* Set azimuth */ Ctrl->A.active = true; j = sscanf(opt->arg, "%lf/%lf", &Ctrl->A.azimuth[0], &Ctrl->A.azimuth[1]); Ctrl->A.two = (j == 2); break; case 'D': /* Find direction of grad|z| */ Ctrl->D.active = true; j = 0; while (opt->arg[j]) { switch (opt->arg[j]) { case 'C': case 'c': Ctrl->D.mode |= 1; break; case 'O': case 'o': Ctrl->D.mode |= 2; break; case 'N': case 'n': Ctrl->D.mode |= 4; break; default: GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -D option: Unrecognized modifier\n"); n_errors++; break; } j++; } break; case 'E': /* Lambertian family radiance */ Ctrl->E.active = true; switch (opt->arg[0]) { case 'p': /* Peucker */ Ctrl->E.mode = 1; break; case 's': /* "simple" Lambertian case */ Ctrl->E.mode = 2; n_errors += GMT_check_condition (GMT, sscanf(&opt->arg[1], "%lf/%lf", &Ctrl->E.azimuth, &Ctrl->E.elevation) != 2, "Syntax error -Es option: Must append azimuth/elevation\n"); break; case 'm': /* Nice algorithm from an old program called manipRaster by Tierry Souriot */ Ctrl->E.mode = 4; j = sscanf(&opt->arg[1], "%lf/%lf", &Ctrl->E.azimuth, &Ctrl->E.elevation); if (j == 0) { /* Use default values */ Ctrl->E.azimuth = 360 - 45; Ctrl->E.elevation = 45; } else if (j == 1) Ctrl->E.elevation = 45; break; default: Ctrl->E.mode = 3; /* "full" Lambertian case */ n_errors += GMT_check_condition (GMT, sscanf(opt->arg, "%lf/%lf", &Ctrl->E.azimuth, &Ctrl->E.elevation) < 2, "Syntax error -E option: Must give at least azimuth and elevation\n"); entry = pos = 0; while (entry < 6 && (GMT_strtok (opt->arg, "/", &pos, ptr))) { switch (entry) { case 2: if (ptr[0] != '=') Ctrl->E.ambient = atof (ptr); break; case 3: if (ptr[0] != '=') Ctrl->E.diffuse = atof (ptr); break; case 4: if (ptr[0] != '=') Ctrl->E.specular = atof (ptr); break; case 5: if (ptr[0] != '=') Ctrl->E.shine = atof (ptr); break; } entry++; } break; } break; case 'G': /* Output grid */ if ((Ctrl->G.active = GMT_check_filearg (GMT, 'G', opt->arg, GMT_OUT))) Ctrl->G.file = strdup (opt->arg); else n_errors++; break; case 'L': /* GMT4 BCs */ if (GMT_compat_check (GMT, 4)) { GMT_Report (API, GMT_MSG_COMPAT, "Warning: Option -L is deprecated; -n+b%s was set instead, use this in the future.\n", opt->arg); strncpy (GMT->common.n.BC, opt->arg, 4U); /* We turn on geographic coordinates if -Lg is given by faking -fg */ /* But since GMT_parse_f_option is private to gmt_init and all it does */ /* in this case are 2 lines bellow we code it here */ if (!strcmp (GMT->common.n.BC, "g")) { GMT_set_geographic (GMT, GMT_IN); GMT_set_geographic (GMT, GMT_OUT); } } else n_errors += GMT_default_error (GMT, opt->option); break; case 'M': /* Geographic data */ if (GMT_compat_check (GMT, 4)) { GMT_Report (API, GMT_MSG_COMPAT, "Warning: Option -M is deprecated; -fg was set instead, use this in the future.\n"); if (!GMT_is_geographic (GMT, GMT_IN)) GMT_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set */ } else n_errors += GMT_default_error (GMT, opt->option); break; case 'N': /* Normalization */ Ctrl->N.active = true; j = 0; if (opt->arg[j]) { if (opt->arg[j] == 't' || opt->arg[j] == 'T') { Ctrl->N.mode = 1; j++; } else if (opt->arg[j] == 'e' || opt->arg[j] == 'E') { Ctrl->N.mode = 2; j++; } n_opt_args = sscanf(&opt->arg[j], "%lf/%lf/%lf", &Ctrl->N.norm, &Ctrl->N.sigma, &Ctrl->N.offset); } break; case 'S': /* Slope grid */ if ((Ctrl->S.active = GMT_check_filearg (GMT, 'S', opt->arg, GMT_OUT))) Ctrl->S.file = strdup (opt->arg); else n_errors++; break; default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } n_errors += GMT_check_condition (GMT, !(Ctrl->A.active || Ctrl->D.active || Ctrl->E.active), "Syntax error: Must specify -A, -D, or -E\n"); n_errors += GMT_check_condition (GMT, Ctrl->S.active && !Ctrl->S.file, "Syntax error -S option: Must specify output file\n"); n_errors += GMT_check_condition (GMT, !Ctrl->G.file && !Ctrl->S.active, "Syntax error -G option: Must specify output file\n"); n_errors += GMT_check_condition (GMT, !Ctrl->In.file, "Syntax error: Must specify input file\n"); n_errors += GMT_check_condition (GMT, Ctrl->N.active && (n_opt_args && Ctrl->N.norm <= 0.0), "Syntax error -N option: Normalization amplitude must be > 0\n"); n_errors += GMT_check_condition (GMT, Ctrl->N.active && (n_opt_args > 1 && Ctrl->N.sigma <= 0.0) , "Syntax error -N option: Sigma must be > 0\n"); n_errors += GMT_check_condition (GMT, Ctrl->E.active && Ctrl->E.mode > 1 && (Ctrl->E.elevation < 0.0 || Ctrl->E.elevation > 90.0), "Syntax error -E option: Use 0-90 degree range for elevation\n"); if (Ctrl->E.active && (Ctrl->A.active || Ctrl->D.active || Ctrl->S.active)) { GMT_Report (API, GMT_MSG_NORMAL, "Warning: -E option overrides -A, -D or -S\n"); Ctrl->A.active = Ctrl->D.active = Ctrl->S.active = false; } return (n_errors ? GMT_PARSE_ERROR : GMT_OK); } #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_grdgradient_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_grdgradient (void *V_API, int mode, void *args) { bool sigma_set = false, offset_set = false, bad, new_grid = false; int p[4], mx, error = 0; unsigned int row, col, n; uint64_t ij, ij0, index, n_used = 0; char format[GMT_BUFSIZ] = {""}, buffer[GMT_GRID_REMARK_LEN160] = {""}; double dx_grid, dy_grid, x_factor, y_factor, dzdx, dzdy, ave_gradient, wesn[4]; double azim, denom, max_gradient = 0.0, min_gradient = 0.0, rpi, lat, output, one; double x_factor2 = 0.0, y_factor2 = 0.0, dzdx2 = 0.0, dzdy2 = 0.0, dzds1, dzds2; double p0 = 0.0, q0 = 0.0, p0q0_cte = 1.0, norm_z, mag, s[3], lim_x, lim_y, lim_z; double k_ads = 0.0, diffuse, spec, r_min = DBL_MAX, r_max = -DBL_MAX, scale; struct GMT_GRID *Surf = NULL, *Slope = NULL, *Out = NULL; struct GRDGRADIENT_CTRL *Ctrl = NULL; struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; struct GMT_OPTION *options = NULL; struct GMTAPI_CTRL *API = GMT_get_API_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */ /*----------------------- Standard module initialization and parsing ----------------------*/ if (API == NULL) return (GMT_NOT_A_SESSION); if (mode == GMT_MODULE_PURPOSE) return (GMT_grdgradient_usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */ options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */ if (!options || options->option == GMT_OPT_USAGE) bailout (GMT_grdgradient_usage (API, GMT_USAGE));/* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_grdgradient_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */ /* Parse the command-line arguments */ GMT = GMT_begin_module (API, THIS_MODULE_LIB, THIS_MODULE_NAME, &GMT_cpy); /* Save current state */ if (GMT_Parse_Common (API, GMT_PROG_OPTIONS, options)) Return (API->error); Ctrl = New_grdgradient_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = GMT_grdgradient_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the grdgradient main code ----------------------------*/ GMT_Report (API, GMT_MSG_VERBOSE, "Processing input grid\n"); GMT_memset (s, 3, double); GMT_set_pad (GMT, 2U); /* Ensure space for BCs in case an API passed pad == 0 */ if (Ctrl->N.active && Ctrl->N.sigma != 0.0) sigma_set = true; if (Ctrl->N.active && Ctrl->N.offset != 0.0) offset_set = true; if (Ctrl->A.active) { /* Get azimuth in 0-360 range */ while (Ctrl->A.azimuth[0] < 0.0) Ctrl->A.azimuth[0] += 360.0; while (Ctrl->A.azimuth[0] > 360.0) Ctrl->A.azimuth[0] -= 360.0; if (Ctrl->A.two) { while (Ctrl->A.azimuth[1] < 0.0) Ctrl->A.azimuth[1] += 360.0; while (Ctrl->A.azimuth[1] > 360.0) Ctrl->A.azimuth[1] -= 360.0; } } if (Ctrl->E.active) { /* Get azimuth in 0-360 range */ while (Ctrl->E.azimuth < 0.0) Ctrl->E.azimuth += 360.0; while (Ctrl->E.azimuth > 360.0) Ctrl->E.azimuth -= 360.0; } if (Ctrl->E.mode == 2) { /* Precalculate constants */ p0 = cosd (90.0 - Ctrl->E.azimuth) * tand (90.0 - Ctrl->E.elevation); q0 = sind (90.0 - Ctrl->E.azimuth) * tand (90.0 - Ctrl->E.elevation); p0q0_cte = sqrt (1.0 + p0*p0 + q0*q0); } else if (Ctrl->E.mode == 3 || Ctrl->E.mode == 4) { /* Precalculate constants */ if (Ctrl->E.mode == 3) { Ctrl->E.elevation = 90 - Ctrl->E.elevation; k_ads = Ctrl->E.ambient + Ctrl->E.diffuse + Ctrl->E.specular; } s[0] = sind (Ctrl->E.azimuth) * cosd (Ctrl->E.elevation); s[1] = cosd (Ctrl->E.azimuth) * cosd (Ctrl->E.elevation); s[2] = sind (Ctrl->E.elevation); } GMT_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Current -R setting, if any */ if ((Surf = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_HEADER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) { Return (API->error); } if (GMT_is_subset (GMT, Surf->header, wesn)) GMT_err_fail (GMT, GMT_adjust_loose_wesn (GMT, wesn, Surf->header), ""); /* Subset requested; make sure wesn matches header spacing */ GMT_grd_init (GMT, Surf->header, options, true); if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_DATA_ONLY, wesn, Ctrl->In.file, Surf) == NULL) { /* Get subset */ Return (API->error); } if (Ctrl->S.active) { /* Want slope grid */ if ((Slope = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_ALLOC, Surf)) == NULL) Return (API->error); } new_grid = GMT_set_outgrid (GMT, Ctrl->In.file, Surf, &Out); /* true if input is a read-only array */ if (GMT_is_geographic (GMT, GMT_IN) && !Ctrl->E.active) { /* Flat-Earth approximation */ dx_grid = GMT->current.proj.DIST_M_PR_DEG * Surf->header->inc[GMT_X] * cosd ((Surf->header->wesn[YHI] + Surf->header->wesn[YLO]) / 2.0); dy_grid = GMT->current.proj.DIST_M_PR_DEG * Surf->header->inc[GMT_Y]; } else { /* Cartesian */ dx_grid = Surf->header->inc[GMT_X]; dy_grid = Surf->header->inc[GMT_Y]; } one = (Ctrl->D.active) ? +1.0 : -1.0; /* With -D we want positive grad direction, not negative as for shading (-A, -E) */ x_factor = one / (2.0 * dx_grid); y_factor = one / (2.0 * dy_grid); if (Ctrl->A.active) { /* Convert azimuths to radians to save multiplication later */ if (Ctrl->A.two) { Ctrl->A.azimuth[1] *= D2R; x_factor2 = x_factor * sin (Ctrl->A.azimuth[1]); y_factor2 = y_factor * cos (Ctrl->A.azimuth[1]); } Ctrl->A.azimuth[0] *= D2R; x_factor *= sin (Ctrl->A.azimuth[0]); y_factor *= cos (Ctrl->A.azimuth[0]); } /* Index offset of 4-star points relative to current node */ mx = Surf->header->mx; /* Need a signed mx for p[3] in line below */ p[0] = 1; p[1] = -1; p[2] = mx; p[3] = -mx; min_gradient = DBL_MAX; max_gradient = -DBL_MAX; ave_gradient = 0.0; if (Ctrl->E.mode == 3) { lim_x = Surf->header->wesn[XHI] - Surf->header->wesn[XLO]; lim_y = Surf->header->wesn[YHI] - Surf->header->wesn[YLO]; lim_z = Surf->header->z_max - Surf->header->z_min; scale = MAX (lim_z, MAX (lim_x, lim_y)); lim_x /= scale; lim_y /= scale; lim_z /= scale; dx_grid /= lim_x; dy_grid /= lim_y; x_factor = -dy_grid / (2.0 * lim_z); y_factor = -dx_grid / (2.0 * lim_z); } for (row = 0, ij0 = 0ULL; row < Surf->header->ny; row++) { /* ij0 is the index in a non-padded grid */ if (GMT_is_geographic (GMT, GMT_IN) && !Ctrl->E.active) { /* Evaluate latitude-dependent factors */ lat = GMT_grd_row_to_y (GMT, row, Surf->header); dx_grid = GMT->current.proj.DIST_M_PR_DEG * Surf->header->inc[GMT_X] * cosd (lat); if (dx_grid > 0.0) x_factor = one / (2.0 * dx_grid); /* Use previous value at the poles */ if (Ctrl->A.active) { if (Ctrl->A.two) x_factor2 = x_factor * sin (Ctrl->A.azimuth[1]); x_factor *= sin (Ctrl->A.azimuth[0]); } } for (col = 0; col < Surf->header->nx; col++, ij0++) { ij = GMT_IJP (Surf->header, row, col); /* Index into padded grid */ for (n = 0, bad = false; !bad && n < 4; n++) if (GMT_is_fnan (Surf->data[ij+p[n]])) bad = true; if (bad) { /* One of star corners = NaN; assign NaN answers and skip to next node */ index = (new_grid) ? ij : ij0; Out->data[index] = GMT->session.f_NaN; if (Ctrl->S.active) Slope->data[ij] = GMT->session.f_NaN; continue; } /* We can now evaluate the central finite differences */ dzdx = (Surf->data[ij+1] - Surf->data[ij-1]) * x_factor; dzdy = (Surf->data[ij-Surf->header->mx] - Surf->data[ij+Surf->header->mx]) * y_factor; if (Ctrl->A.two) { dzdx2 = (Surf->data[ij+1] - Surf->data[ij-1]) * x_factor2; dzdy2 = (Surf->data[ij-Surf->header->mx] - Surf->data[ij+Surf->header->mx]) * y_factor2; } /* Write output to unused NW corner */ if (Ctrl->A.active) { /* Directional derivatives requested */ if (Ctrl->A.two) { dzds1 = dzdx + dzdy; dzds2 = dzdx2 + dzdy2; output = (fabs (dzds1) > fabs (dzds2)) ? dzds1 : dzds2; } else output = dzdx + dzdy; ave_gradient += output; min_gradient = MIN (min_gradient, output); max_gradient = MAX (max_gradient, output); } else if (Ctrl->D.active) { if (dzdx == 0.0 && dzdy == 0.0) /* Flat, so no preferred direction */ azim = GMT->session.d_NaN; else { azim = (Ctrl->D.mode & 1) ? atan2d (dzdy, dzdx) : 90.0 - atan2d (dzdy, dzdx); if (Ctrl->D.mode & 4) azim += 90.0; if (azim < 0.0) azim += 360.0; if (azim >= 360.0) azim -= 360.0; if (Ctrl->D.mode & 2 && azim >= 180) azim -= 180.0; } output = azim; if (Ctrl->S.active) Slope->data[ij] = (float)hypot (dzdx, dzdy); } else { /* Ctrl->E.active */ if (Ctrl->E.mode == 2) output = (1.0 + p0 * dzdx + q0 * dzdy) / (sqrt (1.0 + dzdx * dzdx + dzdy * dzdy) * p0q0_cte); else if (Ctrl->E.mode == 3) { norm_z = dx_grid * dy_grid; mag = d_sqrt (dzdx * dzdx + dzdy * dzdy + norm_z * norm_z); dzdx /= mag; dzdy /= mag; norm_z /= mag; diffuse = MAX (0, s[0] * dzdx + s[1] * dzdy + s[2] * norm_z); spec = specular (dzdx, dzdy, norm_z, s); spec = pow (spec, Ctrl->E.shine); output = (Ctrl->E.ambient + Ctrl->E.diffuse * diffuse + Ctrl->E.specular * spec) / k_ads; } else if (Ctrl->E.mode == 4) output = (dzdy*s[0] + dzdx*s[1] + 2*s[2]) / (sqrt(dzdy * dzdy + dzdx * dzdx + 4)); else /* Peucker method */ output = -0.4285 * (dzdx - dzdy) - 0.0844 * fabs (dzdx + dzdy) + 0.6599; r_min = MIN (r_min, output); r_max = MAX (r_max, output); } index = (new_grid) ? ij : ij0; Out->data[index] = (float)output; n_used++; } } if (!new_grid) { /* We got away with using the input grid by ignoring the pad. Now we must put the pad back in */ GMT_memset (Out->header->pad, 4, int); /* Must set pad to zero first otherwise we cannot add the pad in */ Out->header->mx = Out->header->nx; Out->header->my = Out->header->ny; /* Since there is no pad */ GMT_grd_pad_on (GMT, Out, GMT->current.io.pad); /* Now reinstate the pad */ } if (GMT_is_geographic (GMT, GMT_IN)) { /* Data is geographic */ double sum; /* If the N or S poles are included then we only want a single estimate at these repeating points */ if (Out->header->wesn[YLO] == -90.0 && Out->header->registration == GMT_GRID_NODE_REG) { /* Average all the multiple N pole estimates */ for (col = 0, ij = GMT_IJP (Out->header, 0, 0), sum = 0.0; col < Out->header->nx; col++, ij++) sum += Out->data[ij]; sum /= Out->header->nx; /* Average gradient */ for (col = 0, ij = GMT_IJP (Out->header, 0, 0); col < Out->header->nx; col++, ij++) Out->data[ij] = (float)sum; } if (Out->header->wesn[YLO] == -90.0 && Out->header->registration == GMT_GRID_NODE_REG) { /* Average all the multiple S pole estimates */ for (col = 0, ij = GMT_IJP (Out->header, Out->header->ny - 1, 0), sum = 0.0; col < Out->header->nx; col++, ij++) sum += Out->data[ij]; sum /= Out->header->nx; /* Average gradient */ for (col = 0, ij = GMT_IJP (Out->header, Out->header->ny - 1, 0); col < Out->header->nx; col++, ij++) Out->data[ij] = (float)sum; } } if (Ctrl->E.active) { /* data must be scaled to the [-1,1] interval, but we'll do it into [-.95, .95] to not get too bright */ scale = 1.0 / (r_max - r_min); GMT_grd_loop (GMT, Out, row, col, ij) { if (GMT_is_fnan (Out->data[ij])) continue; Out->data[ij] = (float)((-1.0 + 2.0 * ((Out->data[ij] - r_min) * scale)) * 0.95); } } if (offset_set) ave_gradient = Ctrl->N.offset; else ave_gradient /= n_used; if (Ctrl->A.active) { /* Report some statistics */ if (Ctrl->N.active) { /* Chose normalization */ if (Ctrl->N.mode == 1) { /* atan transformation */ if (sigma_set) denom = 1.0 / Ctrl->N.sigma; else { denom = 0.0; GMT_grd_loop (GMT, Out, row, col, ij) { if (!GMT_is_fnan (Out->data[ij])) denom += pow (Out->data[ij] - ave_gradient, 2.0); } denom = sqrt ((n_used - 1) / denom); Ctrl->N.sigma = 1.0 / denom; } rpi = 2.0 * Ctrl->N.norm / M_PI; GMT_grd_loop (GMT, Out, row, col, ij) { if (!GMT_is_fnan (Out->data[ij])) Out->data[ij] = (float)(rpi * atan ((Out->data[ij] - ave_gradient) * denom)); } Out->header->z_max = rpi * atan ((max_gradient - ave_gradient) * denom); Out->header->z_min = rpi * atan ((min_gradient - ave_gradient) * denom); } else if (Ctrl->N.mode == 2) { /* Exp transformation */ if (!sigma_set) { Ctrl->N.sigma = 0.0; GMT_grd_loop (GMT, Out, row, col, ij) { if (!GMT_is_fnan (Out->data[ij])) Ctrl->N.sigma += fabsf (Out->data[ij]); } Ctrl->N.sigma = M_SQRT2 * Ctrl->N.sigma / n_used; } denom = M_SQRT2 / Ctrl->N.sigma; GMT_grd_loop (GMT, Out, row, col, ij) { if (GMT_is_fnan (Out->data[ij])) continue; if (Out->data[ij] < ave_gradient) { Out->data[ij] = (float)(-Ctrl->N.norm * (1.0 - exp ( (Out->data[ij] - ave_gradient) * denom))); } else { Out->data[ij] = (float)( Ctrl->N.norm * (1.0 - exp (-(Out->data[ij] - ave_gradient) * denom))); } } Out->header->z_max = Ctrl->N.norm * (1.0 - exp (-(max_gradient - ave_gradient) * denom)); Out->header->z_min = -Ctrl->N.norm * (1.0 - exp ( (min_gradient - ave_gradient) * denom)); } else { /* Linear transformation */ if ((max_gradient - ave_gradient) > (ave_gradient - min_gradient)) denom = Ctrl->N.norm / (max_gradient - ave_gradient); else denom = Ctrl->N.norm / (ave_gradient - min_gradient); GMT_grd_loop (GMT, Out, row, col, ij) { if (!GMT_is_fnan (Out->data[ij])) Out->data[ij] = (float)((Out->data[ij] - ave_gradient) * denom); } Out->header->z_max = (max_gradient - ave_gradient) * denom; Out->header->z_min = (min_gradient - ave_gradient) * denom; } } } GMT_set_pad (GMT, API->pad); /* Reset to session default pad before output */ /* Now we write out: */ if (Ctrl->A.active) { if (Ctrl->N.active) strcpy (buffer, "Normalized directional derivative(s)"); else strcpy (buffer, "Directional derivative(s)"); sprintf (format, "\t%s\t%s\t%s\t%s\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out); GMT_Report (API, GMT_MSG_VERBOSE, " Min Mean Max sigma intensities:"); GMT_Report (API, GMT_MSG_VERBOSE, format, min_gradient, ave_gradient, max_gradient, Ctrl->N.sigma); } else { if (Ctrl->E.mode > 1) strcpy (buffer, "Lambertian radiance"); else if (Ctrl->E.mode == 1) strcpy (buffer, "Peucker piecewise linear radiance"); else strcpy (buffer, "Directions of grad (z) [uphill direction]"); } if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out)) Return (API->error); if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_REMARK, buffer, Out)) Return (API->error); if (Ctrl->G.active && GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, Ctrl->G.file, Out) != GMT_OK) { Return (API->error); } if (Ctrl->S.active) { if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Slope)) Return (API->error); if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_REMARK, "Magnitude of grad (z)", Slope)) Return (API->error); if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, Ctrl->S.file, Slope) != GMT_OK) { Return (API->error); } } Return (EXIT_SUCCESS); }