diff --git a/doc/rst/source/grdview.rst b/doc/rst/source/grdview.rst index d3f19cfa48e..45425e56213 100644 --- a/doc/rst/source/grdview.rst +++ b/doc/rst/source/grdview.rst @@ -19,7 +19,7 @@ Synopsis [ |-I|\ [*file*\|\ *intens*\|\ **+a**\ *azimuth*][**+d**][**+m**\ *ambient*][**+n**\ *args*] ] [ |-Jz|\ \|\ **Z**\ *parameters* ] [ |-N|\ [*level*]\ [**+g**\ *fill*] ] -[ |-Q|\ **c**\|\ **i**\|\ **m**\ [**x**\|\ **y**]\|\ **s**\ [**m**]\ [*color*][**+m**] ] +[ |-Q|\ **c**\|\ **i**\|\ **m**\ [**x**\|\ **y**]\|\ **s**\ [**m**]\|\ **g**\ [**m**]\ [*color*][**+m**] ] [ |SYN_OPT-Rz| ] [ |-S|\ *smooth* ] [ |-T|\ [**+o**\ [*pen*]][**+s**] ] @@ -101,7 +101,7 @@ Optional Arguments .. _-Q: -**-Q**\ **c**\|\ **i**\|\ **m**\ [**x**\|\ **y**]\|\ **s**\ [**m**]\ [*color*][**+m**] +**-Q**\ **c**|\ **g**[**m**][**a**]|\ **i**|\ **m**\ [**x**|\ **y**]|\ **s**\ [**m**]\ [*color*][**+m**] Select one of following directives. For any of these choices: - **c** - Image plot, but will make nodes with *z* = NaN transparent, using the color-masking @@ -110,9 +110,15 @@ Optional Arguments - **i** - Image plot. Optionally append the effective dots-per-unit resolution for the rasterization [Default is :term:`GMT_GRAPHICS_DPU`]. - **m** - Mesh plot [Default]. Optionally append *color* for a different mesh paint [white]. - For waterfall plots, append **x** for row or **y** for column profiles). Specify color as for plain **m**. + For waterfall plots, append **x** for row or **y** for column profiles. Specify color as for plain **m**. - **s** - Surface plot. Optionally append **m** to have mesh lines drawn on top of surface. See **-Wm** for setting a specific mesh *pen*. + - **g** - Gouraud-shaded surface with smooth vertex-based color gradients. Optionally append **m** to draw mesh lines + on top of the surface. Append **a** to use alternate diagonal when splitting tiles into triangles. + Gouraud shading produces smoother color transitions than **-Qs** and generates significantly smaller + and faster PostScript files. The **-Qs** option should be used though when precise color transitions + are required at contour levels (to do a `contourf` type figure. See **-Wc**) because the Gouraud-shaded + doesn't cut the tiles along the contour lines as **-Qs** does. A modifier can adjust the color further: diff --git a/src/grdview.c b/src/grdview.c index a03bfb45ace..c4d8977c8e1 100644 --- a/src/grdview.c +++ b/src/grdview.c @@ -102,9 +102,11 @@ struct GRDVIEW_CTRL { bool active, special; bool mask; bool monochrome; + bool gouraud; /* Enable vertex-based color gradients */ bool cpt; int outline; - unsigned int mode; /* GRDVIEW_MESH, GRDVIEW_SURF, GRDVIEW_IMAGE */ + int mode; /* GRDVIEW_MESH, GRDVIEW_SURF, GRDVIEW_IMAGE */ + int diagonal; /* 0=0-2, 1=adaptive */ double dpi; struct GMT_FILL fill; } Q; @@ -290,6 +292,84 @@ GMT_LOCAL void grdview_add_node (bool used[], double x[], double y[], double z[] (*k)++; } +GMT_LOCAL void grdview_paint_gouraud_tile(struct GMT_CTRL *GMT, struct PSL_CTRL *PSL, struct GMT_PALETTE *P, + struct GMT_GRID *I, double *xmesh, double *ymesh, + gmt_grdfloat *Z_vert, uint64_t ij, int *ij_inc, + bool use_intensity, bool monochrome, int diagonal) { + /* Paint a single grid tile using vertex-based Gouraud shading + * xmesh, ymesh: projected 2D coordinates of 4 tile corners [already projected] + * Z_vert: z-values at 4 tile corners + * ij: index of lower-left corner in grid + * ij_inc: offsets to access 4 tile corners [0, 1, 1-mx, -mx] + * use_intensity: apply illumination from grid I + * monochrome: convert to grayscale + * diagonal: 0=0-2, 1=adaptive + */ + double rgb_vert[12], rgb_tri[9]; + double x_tri[3], y_tri[3]; + double intens; + unsigned int k, indices[6]; + struct GMT_PALETTE_HIDDEN *PH; + + /* Get color for each of 4 vertices */ + for (k = 0; k < 4; k++) { + int index = gmt_get_rgb_from_z(GMT, P, Z_vert[k], &rgb_vert[k*3]); + if (k == 0) + PH = gmt_get_C_hidden(P); + + if (use_intensity) { + intens = I->data[ij + ij_inc[k]]; /* Actual vertex intensity */ + gmt_illuminate(GMT, intens, &rgb_vert[k*3]); + } + + if (monochrome) { + double *rgb = &rgb_vert[k*3]; + double gray = gmt_M_yiq (rgb); + rgb[0] = rgb[1] = rgb[2] = gray; + } + } + + /* Determine triangle vertex indices based on diagonal choice */ + if (diagonal == 0) { /* Use 0-2 diagonal */ + indices[0] = 0; indices[1] = 1; indices[2] = 2; /* Triangle 1 */ + indices[3] = 0; indices[4] = 2; indices[5] = 3; /* Triangle 2 */ + } + else { /* Adaptive - choose based on z-variance */ + double var_02 = fabs(Z_vert[0] - Z_vert[2]); + double var_13 = fabs(Z_vert[1] - Z_vert[3]); + if (var_13 < var_02) { /* Use 1-3 diagonal */ + indices[0] = 0; indices[1] = 1; indices[2] = 3; + indices[3] = 1; indices[4] = 2; indices[5] = 3; + } + else { /* Use 0-2 diagonal (default) */ + indices[0] = 0; indices[1] = 1; indices[2] = 2; + indices[3] = 0; indices[4] = 2; indices[5] = 3; + } + } + + /* Render Triangle 1 */ + for (k = 0; k < 3; k++) { + unsigned int idx = indices[k]; + x_tri[k] = xmesh[idx]; + y_tri[k] = ymesh[idx]; + rgb_tri[k*3+0] = rgb_vert[idx*3+0]; + rgb_tri[k*3+1] = rgb_vert[idx*3+1]; + rgb_tri[k*3+2] = rgb_vert[idx*3+2]; + } + PSL_plotgradienttriangle_gouraud(PSL, x_tri, y_tri, rgb_tri); + + /* Render Triangle 2 */ + for (k = 0; k < 3; k++) { + unsigned int idx = indices[3+k]; + x_tri[k] = xmesh[idx]; + y_tri[k] = ymesh[idx]; + rgb_tri[k*3+0] = rgb_vert[idx*3+0]; + rgb_tri[k*3+1] = rgb_vert[idx*3+1]; + rgb_tri[k*3+2] = rgb_vert[idx*3+2]; + } + PSL_plotgradienttriangle_gouraud(PSL, x_tri, y_tri, rgb_tri); +} + GMT_LOCAL void grdview_paint_it (struct GMT_CTRL *GMT, struct PSL_CTRL *PSL, struct GMT_PALETTE *P, double *x, double *y, int n, double z, bool intens, bool monochrome, double intensity, int outline) { int index; double rgb[4]; @@ -457,7 +537,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, -2, "Draw a horizontal plane at z = [minimum grid (or -R) value]. For rectangular projections, append +g " "to paint the facade between the plane and the data perimeter."); GMT_Option (API, "O,P"); - GMT_Usage (API, 1, "\n-Qc|i|m[x|y]|s[][+m]"); + GMT_Usage (API, 1, "\n-Qc|g[m][a]|i|m[x|y]|s[][+m]"); GMT_Usage (API, -2, "Set plot type request via on of the following directives:"); GMT_Usage (API, 3, "c: Transform polygons to raster-image via scanline conversion. Append effective [%lg%c]. " "Set PS Level 3 color-masking for nodes with z = NaN.", @@ -468,6 +548,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { "Alternatively, append x or y for row or column \"waterfall\" profiles.", gmt_putcolor (API->GMT, API->GMT->PSL->init.page_rgb)); GMT_Usage (API, 3, "s: Colored or shaded surface. Optionally, append m to draw mesh-lines on the surface."); + GMT_Usage (API, 3, "g: Gouraud-shaded surface with smooth vertex-based color gradients. Append m to draw mesh-lines. Append a for alternate diagonal."); GMT_Usage (API, -2, "Color may be further adjusted by a modifier:"); GMT_Usage (API, 3, "+m Force a monochrome image using the YIQ transformation."); GMT_Option (API, "R"); @@ -656,6 +737,13 @@ static int parse (struct GMT_CTRL *GMT, struct GRDVIEW_CTRL *Ctrl, struct GMT_OP Ctrl->Q.special = true; Ctrl->Q.cpt = true; /* Will need a CPT */ /* Intentionally fall through - to 'i' */ + case 'g': /* Gouraud-shaded surface */ + Ctrl->Q.mode = GRDVIEW_SURF; + Ctrl->Q.gouraud = true; + Ctrl->Q.cpt = true; /* Will need a CPT */ + if (opt->arg[1] == 'm') Ctrl->Q.outline = 1; + if (opt->arg[1] == 'a') Ctrl->Q.diagonal = 1; + break; case 'i': /* Image with clipmask */ Ctrl->Q.mode = GRDVIEW_IMAGE; if (opt->arg[1] && isdigit ((int)opt->arg[1])) Ctrl->Q.dpi = grdview_get_dpi (&opt->arg[1]); @@ -679,7 +767,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDVIEW_CTRL *Ctrl, struct GMT_OP if (opt->arg[n+1]) { /* Appended / or just */ k = ((opt->arg[n+1] == '/') ? 2 : 1) + n; n_errors += gmt_M_check_condition (GMT, gmt_getfill (GMT, &opt->arg[k], &Ctrl->Q.fill), - "Option -Qm: To give mesh color, use -Qm[x|y]\n"); + "Option -Qm: To give mesh color, use -Qm[x|y]\n"); } break; case 's': /* Color without contours */ @@ -694,7 +782,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDVIEW_CTRL *Ctrl, struct GMT_OP } if (c != NULL && Ctrl->Q.monochrome) c[0] = '+'; /* Restore the chopped off +m */ - else if (gmt_M_compat_check (GMT, 4) && opt->arg[strlen(opt->arg)-1] == 'g') { + else if (gmt_M_compat_check (GMT, 4) && !Ctrl->Q.gouraud && opt->arg[strlen(opt->arg)-1] == 'g') { GMT_Report (API, GMT_MSG_COMPAT, "Option -Q[g] is deprecated; use -Q[+m] in the future.\n"); Ctrl->Q.monochrome = true; } @@ -991,7 +1079,7 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { } if (cpt) gmt_M_str_free (cpt); } - get_contours = (Ctrl->Q.mode == GRDVIEW_MESH && Ctrl->W.contour) || (Ctrl->Q.mode == GRDVIEW_SURF && P && P->n_colors > 1); + get_contours = (Ctrl->Q.mode == GRDVIEW_MESH && Ctrl->W.contour) || (Ctrl->Q.mode == GRDVIEW_SURF && !Ctrl->Q.gouraud && P && P->n_colors > 1) || (Ctrl->Q.gouraud && Ctrl->W.contour); if (Ctrl->G.active) { /* Draping wanted */ if (Ctrl->G.n == 1 && gmt_M_file_is_image (Ctrl->G.file[0])) { @@ -1699,7 +1787,7 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { double *xcont = NULL, *ycont = NULL, *zcont = NULL, *vcont = NULL, X_vert[4], Y_vert[4], saddle_small; gmt_grdfloat Z_vert[4]; - GMT_Report (API, GMT_MSG_INFORMATION, "Place filled surface\n"); + GMT_Report(API, GMT_MSG_INFORMATION, "Place filled surface\n"); /* PW: Bugs fixed in Nov, 2011: Several problems worth remembering: 1) Earlier [2004] we had fixed grdcontour but not grdview in dealing with the current zero contour. Because of gmt_grdfloat precision we cannot take the grid and repeatedly subtract the difference in contour values. @@ -1723,21 +1811,21 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { of the 4 nodes as the other 3 will pull the average into the middle somewhere. */ - xcont = gmt_M_memory (GMT, NULL, max_alloc, double); - ycont = gmt_M_memory (GMT, NULL, max_alloc, double); - zcont = gmt_M_memory (GMT, NULL, max_alloc, double); - vcont = gmt_M_memory (GMT, NULL, max_alloc, double); + xcont = gmt_M_memory(GMT, NULL, max_alloc, double); + ycont = gmt_M_memory(GMT, NULL, max_alloc, double); + zcont = gmt_M_memory(GMT, NULL, max_alloc, double); + vcont = gmt_M_memory(GMT, NULL, max_alloc, double); - PSL_comment (PSL, "Start of filled surface\n"); - if (Ctrl->Q.outline) gmt_setpen (GMT, &Ctrl->W.pen[1]); + PSL_comment(PSL, "Start of filled surface\n"); + if (Ctrl->Q.outline) gmt_setpen(GMT, &Ctrl->W.pen[1]); for (j = start[0]; j != stop[0]; j += inc[0]) { for (i = start[1]; i != stop[1]; i += inc[1]) { if (id[0] == GMT_Y) { y_bottom = yval[j]; x_left = xval[i]; - bin = gmt_M_ij0 (Topo->header, j, i); - ij = gmt_M_ijp (Topo->header, j, i); + bin = gmt_M_ij0(Topo->header, j, i); + ij = gmt_M_ijp(Topo->header, j, i); } else { y_bottom = yval[i]; @@ -1747,14 +1835,14 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { } y_top = y_bottom + Z->header->inc[GMT_Y]; x_right = x_left + Z->header->inc[GMT_X]; - for (k = bad = 0; !bad && k < 4; k++) bad += gmt_M_is_fnan (Topo->data[ij+ij_inc[k]]); + for (k = bad = 0; !bad && k < 4; k++) bad += gmt_M_is_fnan(Topo->data[ij+ij_inc[k]]); if (bad) { if (P->bfn[GMT_NAN].skip || GMT->current.proj.three_D) continue; X_vert[0] = X_vert[3] = x_left; X_vert[1] = X_vert[2] = x_right; Y_vert[0] = Y_vert[1] = y_bottom; Y_vert[2] = Y_vert[3] = y_top; - for (k = 0; k < 4; k++) gmt_geoz_to_xy (GMT, X_vert[k], Y_vert[k], 0.0, &xmesh[k], &ymesh[k]); - grdview_paint_it (GMT, PSL, P, xmesh, ymesh, 4, GMT->session.d_NaN, false, Ctrl->Q.monochrome, 0.0, Ctrl->Q.outline); + for (k = 0; k < 4; k++) gmt_geoz_to_xy(GMT, X_vert[k], Y_vert[k], 0.0, &xmesh[k], &ymesh[k]); + grdview_paint_it(GMT, PSL, P, xmesh, ymesh, 4, GMT->session.d_NaN, false, Ctrl->Q.monochrome, 0.0, Ctrl->Q.outline); continue; } @@ -1767,7 +1855,7 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { this_intensity = Ctrl->I.value; } - PSL_comment (PSL, "Filled surface bin (%d, %d)\n", j, i); + PSL_comment(PSL, "Filled surface bin (%d, %d)\n", j, i); /* Get mesh polygon */ X_vert[0] = X_vert[3] = x_left; X_vert[1] = X_vert[2] = x_right; @@ -1777,7 +1865,6 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { * the contouring stage we need to make the same adjustments below */ for (k = 0; k < 4; k++) Z_vert[k] = Z->data[ij+ij_inc[k]]; /* First a straight copy */ - if (get_contours && binij[bin].first_cont) { /* Contours go through here */ /* Determine if this bin will give us saddle trouble */ @@ -2038,25 +2125,55 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { } else { /* No Contours */ - /* For stability, take the color corresponding to the average value of the four corners */ - z_ave = 0.25 * (Z_vert[0] + Z_vert[1] + Z_vert[2] + Z_vert[3]); + for (k = 0; k < 4; k++) + gmt_geoz_to_xy(GMT, X_vert[k], Y_vert[k], (double)(Topo->data[ij+ij_inc[k]]), &xmesh[k], &ymesh[k]); - /* Now paint the polygon piece */ + if (Ctrl->Q.gouraud) { + /* Gouraud shading - vertex-based colors */ + grdview_paint_gouraud_tile(GMT, PSL, P, Intens, xmesh, ymesh, Z_vert, ij, ij_inc, + Ctrl->I.active, Ctrl->Q.monochrome, Ctrl->Q.diagonal); + /* Draw contour lines if desired */ + pen_set = false; + for (this_cont = start_cont; Ctrl->W.contour && this_cont; this_cont = this_cont->next_cont) { + for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) { + z_val = (Ctrl->G.active) ? gmt_bcr_get_z (GMT, Topo, (double)this_point->x, (double)this_point->y) : this_cont->value; + if (gmt_M_is_dnan (z_val)) continue; - for (k = 0; k < 4; k++) gmt_geoz_to_xy (GMT, X_vert[k], Y_vert[k], (double)(Topo->data[ij+ij_inc[k]]), &xmesh[k], &ymesh[k]); - grdview_paint_it (GMT, PSL, P, xmesh, ymesh, 4, z_ave+small, Ctrl->I.active, Ctrl->Q.monochrome, this_intensity, Ctrl->Q.outline); + gmt_geoz_to_xy (GMT, (double)this_point->x, (double)this_point->y, z_val, &xx[k], &yy[k]); + k++; + } + if (!pen_set) { + gmt_setpen (GMT, &Ctrl->W.pen[0]); + pen_set = true; + } + PSL_plotline (PSL, xx, yy, k, PSL_MOVE|PSL_STROKE); + } + if (pen_set) gmt_setpen (GMT, &Ctrl->W.pen[1]); + if (Ctrl->Q.outline) { + PSL_setfill (PSL, GMT->session.no_rgb, 1); + PSL_plotpolygon (PSL, xmesh, ymesh, 4); + } + } + else { + /* Traditional flat shading */ + /* For stability, take the color corresponding to the average value of the four corners */ + z_ave = 0.25 * (Z_vert[0] + Z_vert[1] + Z_vert[2] + Z_vert[3]); + + /* Now paint the polygon piece */ + grdview_paint_it(GMT, PSL, P, xmesh, ymesh, 4, z_ave+small, Ctrl->I.active, Ctrl->Q.monochrome, this_intensity, Ctrl->Q.outline); + } } } } - gmt_M_free (GMT, xcont); - gmt_M_free (GMT, ycont); - gmt_M_free (GMT, zcont); - gmt_M_free (GMT, vcont); + gmt_M_free(GMT, xcont); + gmt_M_free(GMT, ycont); + gmt_M_free(GMT, zcont); + gmt_M_free(GMT, vcont); } - PSL_setdash (PSL, NULL, 0); + PSL_setdash(PSL, NULL, 0); - if (GMT->current.proj.z_pars[0] == 0.0) gmt_map_clip_off (GMT); + if (GMT->current.proj.z_pars[0] == 0.0) gmt_map_clip_off(GMT); if (Ctrl->N.facade) { /* Cover the two front sides */ PSL_comment (PSL, "Painting the frontal facade\n"); @@ -2105,13 +2222,13 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { } } - gmt_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level); - gmt_map_basemap (GMT); /* Plot basemap last if not 3-D */ + gmt_plane_perspective(GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level); + gmt_map_basemap(GMT); /* Plot basemap last if not 3-D */ if (GMT->current.proj.three_D) - gmt_vertical_axis (GMT, 2); /* Draw foreground axis */ - gmt_plane_perspective (GMT, -1, 0.0); + gmt_vertical_axis(GMT, 2); /* Draw foreground axis */ + gmt_plane_perspective(GMT, -1, 0.0); - gmt_plotend (GMT); + gmt_plotend(GMT); /* Free memory */ @@ -2136,33 +2253,33 @@ EXTERN_MSC int GMT_grdview(void *V_API, int mode, void *args) { gmt_M_free (GMT, binij); } - gmt_change_grdreg (GMT, Topo->header, t_reg); /* Reset registration, if required */ + gmt_change_grdreg(GMT, Topo->header, t_reg); /* Reset registration, if required */ if (use_intensity_grid) { - gmt_change_grdreg (GMT, Intens->header, i_reg); /* Reset registration, if required */ + gmt_change_grdreg(GMT, Intens->header, i_reg); /* Reset registration, if required */ if (saved_data_pointer) { - gmt_M_free (GMT, Intens->data); + gmt_M_free(GMT, Intens->data); Intens->data = saved_data_pointer; } } - gmt_M_free (GMT, xx); - gmt_M_free (GMT, yy); - gmt_M_free (GMT, x); - gmt_M_free (GMT, y); - gmt_M_free (GMT, z); - gmt_M_free (GMT, v); + gmt_M_free(GMT, xx); + gmt_M_free(GMT, yy); + gmt_M_free(GMT, x); + gmt_M_free(GMT, y); + gmt_M_free(GMT, z); + gmt_M_free(GMT, v); if (Ctrl->G.active) { for (k = 0; k < Ctrl->G.n; k++) { - gmt_change_grdreg (GMT, Drape[k]->header, d_reg[k]); /* Reset registration, if required */ + gmt_change_grdreg(GMT, Drape[k]->header, d_reg[k]); /* Reset registration, if required */ } } - if (get_contours && GMT_Destroy_Data (API, &Z) != GMT_NOERROR) { - GMT_Report (API, GMT_MSG_ERROR, "Failed to free Z\n"); + if (get_contours && GMT_Destroy_Data(API, &Z) != GMT_NOERROR) { + GMT_Report(API, GMT_MSG_ERROR, "Failed to free Z\n"); } if (Ctrl->C.active) { - if (Ctrl->C.savecpt && GMT_Write_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, 0, NULL, Ctrl->C.savecpt, P) != GMT_NOERROR) { - GMT_Report (API, GMT_MSG_ERROR, "Failed to save the used CPT in file: %s\n", Ctrl->C.savecpt); + if (Ctrl->C.savecpt && GMT_Write_Data(API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, 0, NULL, Ctrl->C.savecpt, P) != GMT_NOERROR) { + GMT_Report(API, GMT_MSG_ERROR, "Failed to save the used CPT in file: %s\n", Ctrl->C.savecpt); } - if (GMT_Destroy_Data (API, &P) != GMT_NOERROR) {Return (API->error);} + if (GMT_Destroy_Data(API, &P) != GMT_NOERROR) {Return (API->error);} } Return (GMT_NOERROR); diff --git a/src/postscriptlight.c b/src/postscriptlight.c index 7bf16adaa42..a1df1ada0e0 100644 --- a/src/postscriptlight.c +++ b/src/postscriptlight.c @@ -4981,8 +4981,9 @@ int PSL_plotgradienttriangle_gouraud(struct PSL_CTRL *PSL, double *x, double *y, * x, y: arrays of 3 coordinates (in inches, will be converted to points) * rgb: array of 9 values [r1 g1 b1 r2 g2 b2 r3 g3 b3], each 0.0-1.0 * This uses PostScript Level 2/3 shading for perfectly smooth gradients - * - * Writen by Claude.ai + * + * Optimized: /GT procedure with proper variable handling + * Written by Claude.ai */ int i; @@ -4994,22 +4995,25 @@ int PSL_plotgradienttriangle_gouraud(struct PSL_CTRL *PSL, double *x, double *y, y_pts[i] = y[i] * PSL->internal.dpu; } + /* Define Gouraud Triangle procedure once */ + static int gt_defined = 0; + if (!gt_defined) { + PSL_command(PSL, "/GT {\n"); + PSL_command(PSL, " /data exch def /ymax exch def /ymin exch def /xmax exch def /xmin exch def\n"); + PSL_command(PSL, " << /ShadingType 4 /ColorSpace /DeviceRGB\n"); + PSL_command(PSL, " /BitsPerCoordinate 16 /BitsPerComponent 8 /BitsPerFlag 8\n"); + PSL_command(PSL, " /Decode [xmin xmax ymin ymax 0 1 0 1 0 1]\n"); + PSL_command(PSL, " /DataSource data >>\n"); + PSL_command(PSL, " shfill\n"); + PSL_command(PSL, "} def\n"); + gt_defined = 1; + } + PSL_comment(PSL, "Begin Gouraud shaded triangle\n"); PSL_command(PSL, "gsave\n"); - /* Set clipping path to triangle to constrain shfill and establish bounding box */ - PSL_command(PSL, "N %.4f %.4f M %.4f %.4f L %.4f %.4f L closepath clip", + PSL_command(PSL, "N %.4f %.4f M %.4f %.4f L %.4f %.4f L closepath clip\n", x_pts[0], y_pts[0], x_pts[1], y_pts[1], x_pts[2], y_pts[2]); - - /* Create Type 4 (Free-Form Gouraud) shading dictionary */ - PSL_command(PSL, "<<\n"); - PSL_command(PSL, " /ShadingType 4\n"); - PSL_command(PSL, " /ColorSpace /DeviceRGB\n"); - PSL_command(PSL, " /BitsPerCoordinate 16\n"); - PSL_command(PSL, " /BitsPerComponent 8\n"); - PSL_command(PSL, " /BitsPerFlag 8\n"); - - /* Decode arrays map from integer coordinates to actual coordinate ranges */ /* Find bounding box */ double xmin = x_pts[0], xmax = x_pts[0], ymin = y_pts[0], ymax = y_pts[0]; for (i = 1; i < 3; i++) { @@ -5019,25 +5023,19 @@ int PSL_plotgradienttriangle_gouraud(struct PSL_CTRL *PSL, double *x, double *y, if (y_pts[i] > ymax) ymax = y_pts[i]; } - PSL_command(PSL, " /Decode [%.4f %.4f %.4f %.4f 0 1 0 1 0 1]\n", xmin, xmax, ymin, ymax); - PSL_command(PSL, " /DataSource <\n"); - - /* Write triangle data as hex string */ - /* Flag 0 = start new triangle, then coordinates (x,y) scaled to 0-65535, then RGB scaled to 0-255 */ + /* Push parameters onto stack: xmin xmax ymin ymax */ + PSL_command(PSL, "%.4f %.4f %.4f %.4f <", xmin, xmax, ymin, ymax); for (i = 0; i < 3; i++) { - unsigned int flag = 0; /* Start new triangle with vertex 0 */ + unsigned int flag = 0; unsigned int x_scaled = (unsigned int)((x_pts[i] - xmin) / (xmax - xmin) * 65535.0 + 0.5); unsigned int y_scaled = (unsigned int)((y_pts[i] - ymin) / (ymax - ymin) * 65535.0 + 0.5); unsigned int r = (unsigned int)(rgb[i*3 + 0] * 255.0 + 0.5); unsigned int g = (unsigned int)(rgb[i*3 + 1] * 255.0 + 0.5); unsigned int b = (unsigned int)(rgb[i*3 + 2] * 255.0 + 0.5); - PSL_command(PSL, " %02X %04X %04X %02X%02X%02X\n", flag, x_scaled, y_scaled, r, g, b); + PSL_command(PSL, " %02X%04X%04X%02X%02X%02X", flag, x_scaled, y_scaled, r, g, b); } - - PSL_command(PSL, " >\n"); - PSL_command(PSL, ">>\n"); - PSL_command(PSL, "shfill\n"); + PSL_command(PSL, "> GT\n"); PSL_command(PSL, "grestore\n"); PSL_comment(PSL, "End Gouraud shaded triangle\n"); diff --git a/test/grdview/gouraud.sh b/test/grdview/gouraud.sh new file mode 100644 index 00000000000..7451a2ac865 --- /dev/null +++ b/test/grdview/gouraud.sh @@ -0,0 +1,27 @@ +#!/usr/bin/env bash +# +# Test Gouraud shading for grdview surfaces with -Qg +# Creates comparison plots: flat shading vs Gouraud shading + +ps=gouraud.ps + +# Create a simple synthetic grid (cone/peak) +gmt grdmath -R-5/5/-5/5 -I0.2 X Y HYPOT NEG 3 POW 10 MUL = peak.nc + +# Create CPT +gmt makecpt -Chot -T-3536/0 > t.cpt + +# Row 1: Flat shading (traditional) +gmt grdview peak.nc -R-5/5/-5/5/-3536/0 -JX8c -JZ4c -Qs -Ct.cpt -p135/30 \ + -Bxafg -Byafg -Bzafg -BWSneZ+t"Flat Shading (-Qs)" -P -K -Xc -Y18c > $ps + +# Row 2: Gouraud shading (default diagonal) +gmt grdview peak.nc -R-5/5/-5/5/-3536/0 -JX8c -JZ4c -Qg -Ct.cpt -p135/30 \ + -Bxafg -Byafg -Bzafg -BWSneZ+t"Gouraud (-Qg)" -O -K -Y-10c >> $ps + +# Row 3: Gouraud with mesh +gmt grdview peak.nc -R-5/5/-5/5/-3536/0 -JX8c -JZ4c -Qgm -Ct.cpt -p135/30 \ + -Bxafg -Byafg -Bzafg -BWSneZ+t"Gouraud with mesh (-Qgm)" -O -Y-10c >> $ps + +# Cleanup +rm -f peak.nc t.cpt