diff --git a/doc/rst/source/grdimage_common.rst_ b/doc/rst/source/grdimage_common.rst_ index 38a7228623c..1b98c4a80c6 100644 --- a/doc/rst/source/grdimage_common.rst_ +++ b/doc/rst/source/grdimage_common.rst_ @@ -4,7 +4,7 @@ Description ----------- **grdimage** reads a 2-D grid file and produces a gray-shaded (or -colored) map by by building a rectangular image and assigning pixels +colored) map by building a rectangular image and assigning pixels a gray-shade (or color) based on the z-value and the CPT file. Optionally, illumination may be added by providing a file with intensities in the (-1,+1) range or instructions to derive intensities diff --git a/src/gmt_grd.h b/src/gmt_grd.h index 2a667c8735a..8ce8081989b 100644 --- a/src/gmt_grd.h +++ b/src/gmt_grd.h @@ -150,7 +150,7 @@ enum gmt_enum_wesnids { /*! IJ macro using h but treats the entire grid with pad as no-pad grid, i.e. using mx as width */ #define gmt_M_ij(h,row,col) ((uint64_t)(((int64_t)(row))*((int64_t)h->mx)+(int64_t)(col))) /*! IJPGI macro using h and the pad info that works for either grids (n_bands = 1) or images (n_bands = 1,3,4) */ -#define gmt_M_ijpgi(h,row,col) ((uint64_t)(((int64_t)(row)+(int64_t)h->pad[YHI])*((int64_t)h->mx*(int64_t)h->n_bands)+(int64_t)(col)+(int64_t)h->pad[XLO]*(int64_t)h->n_bands)) +#define gmt_M_ijpgi(h,row,col) ((uint64_t)(((int64_t)(row)+(int64_t)h->pad[YHI])*((int64_t)h->mx*(int64_t)h->n_bands)+((int64_t)(col)+(int64_t)h->pad[XLO])*(int64_t)h->n_bands)) /*! Obtain row and col from index */ #define gmt_M_col(h,ij) (((ij) % h->mx) - h->pad[XLO]) diff --git a/src/grdimage.c b/src/grdimage.c index d239f64f2cb..d833972ea45 100644 --- a/src/grdimage.c +++ b/src/grdimage.c @@ -41,6 +41,8 @@ static char *gdal_ext[N_IMG_EXTENSIONS] = {"tiff", "tif", "gif", "png", "jpg", "bmp"}; #endif +#define GRDIMAGE_NAN_INDEX (GMT_NAN - 3) + /* Control structure for grdimage */ struct GRDIMAGE_CTRL { @@ -100,6 +102,21 @@ struct GRDIMAGE_CTRL { } W; }; +struct GRDIMAGE_CONF { + /* Configuration structure for things to pass around to sub-functions */ + unsigned int colormask_offset; /* Either 0 or 3 depending on -Q */ + bool int_mode; /* true if we are to apply illumination */ + unsigned int *actual_row; /* Array with actual rows as a function of pseudo row */ + unsigned int *actual_col; /* Array of actual columns as a function of pseudo col */ + int n_columns, n_rows; /* Signed dimensions for use in OpenMP loops */ + uint64_t nm; /* Number of pixels in the image */ + struct GMT_PALETTE *P; /* Pointer to the active palette [NULL if image] */ + struct GMT_GRID_HEADER *H; /* Pointer to the active header */ + struct GMT_GRID *Grid; /* Pointer to the active grid [NULL if image] */ + struct GMT_GRID *Intens; /* Pointer to the active intensity grid [NULL if no intensity] */ + struct GMT_IMAGE *Image; /* Pointer to the active image [NUYLL if grid] */ +}; + static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ struct GRDIMAGE_CTRL *C; @@ -399,7 +416,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GMT_O case 'Q': /* PS3 colormasking */ Ctrl->Q.active = true; break; - case 'W': /* Warn if no image */ + case 'W': /* Warn if no image, usually when called from grd2kml */ Ctrl->W.active = true; break; @@ -591,6 +608,425 @@ GMT_LOCAL void grdimage_set_proj_limits (struct GMT_CTRL *GMT, struct GMT_GRID_H } } +GMT_LOCAL void grdimage_blackwhite_PS_image (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, unsigned char *image, unsigned int n_columns, unsigned int n_rows, double x_LL, double y_LL, double dx, double dy) { + /* This function takes an 8-bit grayshade image that we know is only white (255) or black (0) and converts it + * to a 1-bit black/white image suitable for PSL to plot using PostScripts image operator for 1-bit images. + * Because all projections and scalings have already taken place, this is a simple scanline operation. */ + unsigned char *bitimage = NULL; + unsigned int row, col, nx8, shift, b_or_w, nx_pixels; + uint64_t imsize, k, k8, byte; + double x_side, y_side = n_rows * dy; + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Converting 8-bit image to 1-bit B/W image then plotting it\n"); + + nx8 = irint (ceil (n_columns / 8.0)); /* Image width must be a multiple of 8 bits, so we round up */ + nx_pixels = nx8 * 8; /* The row length in bits after the rounding up */ + imsize = gmt_M_get_nm (GMT, nx8, n_rows); + bitimage = gmt_M_memory (GMT, NULL, imsize, unsigned char); /* Memory to hold the 1-bit image */ + + /* Reprocess the byte image. Here there are no worries about direction of rows, cols since that was dealt with during color assignment */ + + for (row = 0, k = k8 = 0; row < n_rows; row++) { /* Process each scan line */ + shift = 0; byte = 0; + for (col = 0; col < n_columns; col++, k++) { /* Visit each byte in the original gray shade image */ + b_or_w = (image[k] == 255); /* Let white == 1, black == 0 */ + byte |= b_or_w; /* Add in the current bit (0 or 1) */ + shift++; /* Position us for the next bit */ + if (shift == 8) { /* Did all 8, time to dump out another byte ["another one bytes the dust", he, he,...] */ + bitimage[k8++] = (unsigned char) byte; /* Place the last 8 bits in output array... */ + byte = shift = 0; /* ...and start over for next sequence of 8 bit nodes */ + } + else /* Move the bits we have so far 1 step to the left */ + byte <<= 1; + } + if (shift) { /* Set the remaining bits in this bit to white; this can apply to the last 1-7 nodes on each row */ + byte |= 1; + shift++; + while (shift < 8) { + byte <<= 1; + byte |= 1; + shift++; + } + bitimage[k8++] = (unsigned char) byte; /* Copy final byte from this row into the image */ + } + } + + x_side = nx_pixels * dx; /* Since the image may be 1-7 bits wider than what we need we may enlarge it by a tiny amount */ + PSL_plotbitimage (GMT->PSL, x_LL, y_LL, x_side, y_side, PSL_BL, bitimage, nx_pixels, n_rows, Ctrl->G.rgb[GMT_FGD], Ctrl->G.rgb[GMT_BGD]); + gmt_M_free (GMT, bitimage); /* Done with the B/W image array */ +} + +GMT_LOCAL void grdimage_grd_gray_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) grid, 2) grayscale, 3) no intensity */ + int srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + gmt_M_unused (Ctrl); + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> gray image with no illumination.\n"); +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,scol,node,rgb) shared(GMT,Conf,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + (void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + image[byte++] = gmt_M_u255 (rgb[0]); /* Color table only has grays, just use r since r = g = b here */ + } + } +} + +GMT_LOCAL void grdimage_grd_gray_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) grid, 2) grayscale, 3) with intensity */ + int srow, scol, index; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> gray image with illumination.\n"); +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,scol,node,index,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + if (index != GRDIMAGE_NAN_INDEX) { /* Add in the effect of illumination */ + if (Conf->int_mode) /* Intensity value comes from a co-registered grid */ + gmt_illuminate (GMT, Conf->Intens->data[node], rgb); + else /* A constant (ambient) intensity was given via -I */ + gmt_illuminate (GMT, Ctrl->I.value, rgb); + } + image[byte++] = gmt_M_u255 (rgb[0]); /* Color table only has grays, just use r since r = g = b here */ + } + } +} + +GMT_LOCAL void grdimage_grd_c2s_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) grid, 2) color -> gray via YIQ, 3) no intensity */ + int srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + gmt_M_unused (Ctrl); + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n"); +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,scol,node,rgb) shared(GMT,Conf,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + (void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb)); + } + } +} + +GMT_LOCAL void grdimage_grd_c2s_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) grid, 2) color -> gray via YIQ, 3) no intensity */ + int srow, scol, index; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n"); +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,scol,node,index,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */ + if (Conf->int_mode) /* Intensity value comes from the grid */ + gmt_illuminate (GMT, Conf->Intens->data[node], rgb); + else /* A constant (ambient) intensity was given via -I */ + gmt_illuminate (GMT, Ctrl->I.value, rgb); + } + image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb)); + } + } +} + +GMT_LOCAL void grdimage_grd_color_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */ + int srow, scol, k; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + gmt_M_unused (Ctrl); + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n"); +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,k,scol,node,rgb) shared(GMT,Conf,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = Conf->colormask_offset + 3 * srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + (void)gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]); + } + } +} + +GMT_LOCAL void grdimage_grd_color_no_intensity_CM (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image, unsigned char *rgb_used) { + /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */ + unsigned char i_rgb[3]; + int srow, scol, k, index; + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + gmt_M_unused (Ctrl); + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n"); + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = Conf->colormask_offset + 3 * srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + for (k = 0; k < 3; k++) image[byte++] = i_rgb[k] = gmt_M_u255 (rgb[k]); + if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */ + index = (i_rgb[0]*256 + i_rgb[1])*256 + i_rgb[2]; /* The index into the cube for the selected NaN color */ + rgb_used[index] = true; + } + } + } +} + +GMT_LOCAL void grdimage_grd_color_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */ + int srow, scol, index, k; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n"); +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,k,scol,node,index,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = Conf->colormask_offset + 3 * srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + if (index != GRDIMAGE_NAN_INDEX) { /* Deal with illumination */ + if (Conf->int_mode) /* Intensity value comes from the grid */ + gmt_illuminate (GMT, Conf->Intens->data[node], rgb); + else /* A constant (ambient) intensity was given via -I */ + gmt_illuminate (GMT, Ctrl->I.value, rgb); + } + for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]); + } + } +} + +GMT_LOCAL void grdimage_grd_color_with_intensity_CM (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image, unsigned char *rgb_used) { + /* Function that fills out the image in the special case of 1) grid, 2) color, 3) no intensity */ + unsigned char i_rgb[3], n_rgb[3]; + int srow, scol, index, k; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Basic z(x,y) -> color image with no illumination.\n"); + /* Determine NaN rgb 0-255 triple ones */ + (void)gmt_get_rgb_from_z (GMT, Conf->P, GMT->session.d_NaN, rgb); + for (k = 0; k < 3; k++) n_rgb[k] = gmt_M_u255 (rgb[k]); + + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines */ + byte = Conf->colormask_offset + 3 * srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Current grid node */ + index = gmt_get_rgb_from_z (GMT, Conf->P, Conf->Grid->data[node], rgb); + if (index == GRDIMAGE_NAN_INDEX) { /* Nan color */ + for (k = 0; k < 3; k++) image[byte++] = n_rgb[k]; + } + else { /* Deal with illumination for non-NaN nodes */ + if (Conf->int_mode) /* Intensity value comes from the grid */ + gmt_illuminate (GMT, Conf->Intens->data[node], rgb); + else /* A constant (ambient) intensity was given via -I */ + gmt_illuminate (GMT, Ctrl->I.value, rgb); + for (k = 0; k < 3; k++) i_rgb[k] = image[byte++] = gmt_M_u255 (rgb[k]); + index = (i_rgb[0]*256 + i_rgb[1])*256 + i_rgb[2]; /* The index into the cube for the selected NaN color */ + rgb_used[index] = true; + } + } + } +} + +GMT_LOCAL void grdimage_img_set_transparency (struct GMT_CTRL *GMT, unsigned char pix4, double *rgb) { + /* JL: Here we assume background color is white, hence t * 1. + But what would it take to have a user selected background color? */ + double o, t; /* o - opacity, t = transparency */ + o = pix4 / 255.0; t = 1 - o; + rgb[0] = o * rgb[0] + t * GMT->current.map.frame.fill[GMT_Z].rgb[0]; + rgb[1] = o * rgb[1] + t * GMT->current.map.frame.fill[GMT_Z].rgb[1]; + rgb[2] = o * rgb[2] + t * GMT->current.map.frame.fill[GMT_Z].rgb[2]; +} + +GMT_LOCAL void grdimage_img_gray_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) image, 2) gray, 3) with intensity */ + int srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */ + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,scol,node,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Start of current pixel node */ + rgb[0] = rgb[1] = rgb[2] = gmt_M_is255 (Conf->Image->data[node++]); + if (Conf->int_mode) { /* Intensity value comes from the grid, so update node */ + node = gmt_M_ijp (Conf->Intens->header, Conf->actual_row[srow], Conf->actual_col[scol]); + gmt_illuminate (GMT, Conf->Intens->data[node], rgb); /* Apply illumination to this color */ + } + else /* A constant (ambient) intensity was given via -I */ + gmt_illuminate (GMT, Ctrl->I.value, rgb); /* Apply constant illumination to this color */ + image[byte++] = gmt_M_u255 (rgb[0]); + } + } +} + +GMT_LOCAL void grdimage_img_gray_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) image, 2) gray, 3) no intensity */ + int srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */ + gmt_M_unused (Ctrl); + uint64_t byte, kk, node; + gmt_M_unused (GMT); + +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,scol,node) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol]; /* Start of current pixel node */ + image[byte++] = Conf->Image->data[node]; + } + } +} + +GMT_LOCAL void grdimage_img_c2s_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) image, 2) color -> gray via YIQ, 3) with intensity */ + bool transparency = (Conf->Image->header->n_bands == 4); + int srow, scol, k; /* Due to OPENMP on Windows requiring signed int loop variables */ + unsigned n_bands = Conf->Image->header->n_bands; + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,k,scol,node,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol] * n_bands; /* Start of current pixel node */ + for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node++]); + if (transparency && Conf->Image->data[node] < 255) /* Dealing with an image with transparency values less than 255 */ + grdimage_img_set_transparency (GMT, Conf->Image->data[node], rgb); + if (Conf->int_mode) { /* Intensity value comes from the grid, so update node */ + node = gmt_M_ijp (Conf->Intens->header, Conf->actual_row[srow], Conf->actual_col[scol]); + gmt_illuminate (GMT, Conf->Intens->data[node], rgb); /* Apply illumination to this color */ + } + else /* A constant (ambient) intensity was given via -I */ + gmt_illuminate (GMT, Ctrl->I.value, rgb); /* Apply constant illumination to this color */ + image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb)); + } + } +} + +GMT_LOCAL void grdimage_img_c2s_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) image, 2) color -> gray via YIQ, 3) with intensity */ + bool transparency = (Conf->Image->header->n_bands == 4); + int srow, scol, k; /* Due to OPENMP on Windows requiring signed int loop variables */ + unsigned n_bands = Conf->Image->header->n_bands; + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + gmt_M_unused (Ctrl); + +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,k,scol,node,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */ + byte = srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol] * n_bands; /* Start of current pixel node */ + for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node++]); + if (transparency && Conf->Image->data[node] < 255) /* Dealing with an image with transparency values less than 255 */ + grdimage_img_set_transparency (GMT, Conf->Image->data[node], rgb); + image[byte++] = gmt_M_u255 (gmt_M_yiq (rgb)); + } + } +} + +GMT_LOCAL void grdimage_img_color_no_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) image, 2) color, 3) with intensity */ + bool transparency = (Conf->Image->header->n_bands == 4); + int srow, scol, k; /* Due to OPENMP on Windows requiring signed int loop variables */ + unsigned n_bands = Conf->Image->header->n_bands; + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + gmt_M_unused (Ctrl); + +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,k,scol,node,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */ + byte = Conf->colormask_offset + 3 * srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol] * n_bands; /* Start of current pixel node */ + for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node++]); + if (transparency && Conf->Image->data[node] < 255) /* Dealing with an image with transparency values less than 255 */ + grdimage_img_set_transparency (GMT, Conf->Image->data[node], rgb); + for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]); /* Scale up to integer 0-255 range */ + } + } +} + +GMT_LOCAL void grdimage_img_color_with_intensity (struct GMT_CTRL *GMT, struct GRDIMAGE_CTRL *Ctrl, struct GRDIMAGE_CONF *Conf, unsigned char *image) { + /* Function that fills out the image in the special case of 1) image, 2) color, 3) with intensity */ + bool transparency = (Conf->Image->header->n_bands == 4); + int srow, scol, k; /* Due to OPENMP on Windows requiring signed int loop variables */ + unsigned n_bands = Conf->Image->header->n_bands; + uint64_t byte, kk, node; + double rgb[4] = {0.0, 0.0, 0.0, 0.0}; + +#ifdef _OPENMP +#pragma omp parallel for private(srow,byte,kk,k,scol,node,rgb) shared(GMT,Conf,Ctrl,image) +#endif + for (srow = 0; srow < Conf->n_rows; srow++) { /* March along scanlines in the output bitimage */ + byte = Conf->colormask_offset + 3 * srow * Conf->n_columns; + kk = gmt_M_ijpgi (Conf->H, Conf->actual_row[srow], 0); /* Start pixel of this row */ + for (scol = 0; scol < Conf->n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ + node = kk + Conf->actual_col[scol] * n_bands; /* Start of current pixel node */ + for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Conf->Image->data[node++]); + if (transparency && Conf->Image->data[node] < 255) /* Dealing with an image with transparency values less than 255 */ + grdimage_img_set_transparency (GMT, Conf->Image->data[node], rgb); + if (Conf->int_mode) { /* Intensity value comes from the grid, so update node */ + node = gmt_M_ijp (Conf->Intens->header, Conf->actual_row[srow], Conf->actual_col[scol]); + gmt_illuminate (GMT, Conf->Intens->data[node], rgb); /* Apply illumination to this color */ + } + else /* A constant (ambient) intensity was given via -I */ + gmt_illuminate (GMT, Ctrl->I.value, rgb); /* Apply constant illumination to this color */ + for (k = 0; k < 3; k++) image[byte++] = gmt_M_u255 (rgb[k]); /* Scale up to integer 0-255 range */ + } + } +} + #define bailout(code) {gmt_M_free_options (mode); return (code);} #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);} @@ -604,17 +1040,15 @@ GMT_LOCAL void grdimage_set_proj_limits (struct GMT_CTRL *GMT, struct GMT_GRID_H EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { bool done, need_to_project, normal_x, normal_y, resampled = false, gray_only = false; bool nothing_inside = false, use_intensity_grid = false, got_data_tiles = false, rgb_cube_scan; - bool has_content = false, mem_G = false, mem_I = false, mem_D = false, got_z_grid = true; - unsigned int n_columns = 0, n_rows = 0, grid_registration = GMT_GRID_NODE_REG, intensity_mode; - unsigned int colormask_offset = 0, try, row, col, mixed = 0, *actual_row = NULL, *actual_col = NULL; - uint64_t node_RGBA = 0; /* uint64_t for the RGB(A) image array. */ - uint64_t node, k, kk, byte, step, dim[GMT_DIM_SIZE] = {0, 0, 3, 0}; - int index = 0, ks, error = 0, ret_val = GMT_NOERROR, ftype = GMT_NOTSET; + bool has_content, mem_G = false, mem_I = false, mem_D = false, got_z_grid = true; + unsigned int grid_registration = GMT_GRID_NODE_REG, try, row, col, mixed = 0; + uint64_t node, k, kk, dim[GMT_DIM_SIZE] = {0, 0, 3, 0}; + int error = 0, ret_val = GMT_NOERROR, ftype = GMT_NOTSET; char *img_ProjectionRefPROJ4 = NULL, *way[2] = {"via GDAL", "directly"}, cmd[GMT_LEN256] = {""}, data_grd[GMT_VF_LEN] = {""}; - unsigned char *bitimage_8 = NULL, *bitimage_24 = NULL, *rgb_used = NULL, i_rgb[3]; + unsigned char *bitimage_8 = NULL, *bitimage_24 = NULL, *rgb_used = NULL; - double dx, dy, x_side, y_side, x0 = 0.0, y0 = 0.0, rgb[4] = {0.0, 0.0, 0.0, 0.0}; + double dx, dy, x_side, y_side, x0 = 0.0, y0 = 0.0; double img_wesn[4], img_inc[2] = {1.0, 1.0}; /* Image increments & min/max for writing images or external interfaces */ double *NaN_rgb = NULL, red[4] = {1.0, 0.0, 0.0, 0.0}, black[4] = {0.0, 0.0, 0.0, 0.0}, wesn[4] = {0.0, 0.0, 0.0, 0.0}; double *Ix = NULL, *Iy = NULL; /* Pointers to hold on to read-only x/y image arrays */ @@ -633,6 +1067,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { struct GMT_IMAGE *I = NULL, *Img_proj = NULL; /* A GMT image datatype, if GDAL is used */ struct GMT_IMAGE *Out = NULL; /* A GMT image datatype, if external interface is used with -A */ struct GMT_GRID *G2 = NULL; + struct GRDIMAGE_CONF *Conf = NULL; /*----------------------- Standard module initialization and parsing ----------------------*/ @@ -651,6 +1086,11 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { /*---------------------------- This is the grdimage main code ----------------------------*/ + if ((Conf = gmt_M_memory (GMT, NULL, 1, struct GRDIMAGE_CONF)) == NULL) { + GMT_Report (API, GMT_MSG_ERROR, "Failed to allocate Conf structure\n"); + Return (GMT_MEMORY_ERROR); + } + gmt_grd_set_datapadding (GMT, true); /* Turn on gridpadding when reading a subset */ use_intensity_grid = (Ctrl->I.active && !Ctrl->I.constant); /* We want to use an intensity grid */ @@ -838,15 +1278,13 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { gmt_plane_perspective (GMT, -1, 0.0); gmt_plotend (GMT); } - Return (GMT_NOERROR); + if (Ctrl->W.active) ret_val = GMT_IMAGE_NO_DATA; /* Flag that output image has no data - needed by grd2kml only so far */ + Return (ret_val); } /* Here the grid/image is inside the plot domain. The same must be true of any * auto-derived intensities we may create below */ - if (Ctrl->W.active && got_z_grid && gmt_M_is_dnan (header_work->z_min)) - ret_val = GMT_IMAGE_NO_DATA; /* Flag that our output image has no information - this is needed by grd2kml only so far */ - if (Ctrl->I.derive) { /* Auto-create intensity grid from data grid using the now determined data region */ bool got_int4_grid = false; char int_grd[GMT_VF_LEN] = {""}, int4_grd[GMT_VF_LEN] = {""}; @@ -891,12 +1329,12 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { if (got_z_grid) { /* Get grid dimensions */ double *region = (gmt_file_is_tiled_list (API, Ctrl->In.file, NULL, NULL, NULL)) ? API->tile_wesn : wesn; /* Make sure we get correct dimensions if tiled grids are used */ - n_columns = gmt_M_get_n (GMT, region[XLO], region[XHI], Grid_orig->header->inc[GMT_X], Grid_orig->header->registration); - n_rows = gmt_M_get_n (GMT, region[YLO], region[YHI], Grid_orig->header->inc[GMT_Y], Grid_orig->header->registration); + Conf->n_columns = gmt_M_get_n (GMT, region[XLO], region[XHI], Grid_orig->header->inc[GMT_X], Grid_orig->header->registration); + Conf->n_rows = gmt_M_get_n (GMT, region[YLO], region[YHI], Grid_orig->header->inc[GMT_Y], Grid_orig->header->registration); } else { /* Trust the image info from GDAL to make it more stable against pixel vs grid registration troubles */ - n_columns = I->header->n_columns; - n_rows = I->header->n_rows; + Conf->n_columns = I->header->n_columns; + Conf->n_rows = I->header->n_rows; } if (!Ctrl->A.active) { /* Initialize the PostScript plot */ @@ -942,7 +1380,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { /* Create a virtual file to hold the resampled grid - out_string then holds the name of this output "file" */ GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, out_string); /* Create the command to do the resampling via the grdsample module */ - sprintf (cmd, "%s -G%s -I%d+/%d+ --GMT_HISTORY=readonly", in_string, out_string, n_columns, n_rows); + sprintf (cmd, "%s -G%s -I%d+/%d+ --GMT_HISTORY=readonly", in_string, out_string, Conf->n_columns, Conf->n_rows); GMT_Report (API, GMT_MSG_INFORMATION, "Calling grdsample with args %s\n", cmd); if (GMT_Call_Module (GMT->parent, "grdsample", GMT_MODULE_CMD, cmd) != GMT_NOERROR) /* Do the resampling */ return (API->error); @@ -953,7 +1391,6 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { } Intens_orig = G2; /* ...and point to the resampled intensity grid */ } - } if (got_z_grid && (gmt_whole_earth (GMT, Grid_orig->header->wesn, wesn) == 1)) @@ -964,8 +1401,8 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { double inc[2] = {0.0, 0.0}; if (Ctrl->E.dpi == 0) { /* Use input # of nodes as # of projected nodes */ - nx_proj = n_columns; - ny_proj = n_rows; + nx_proj = Conf->n_columns; + ny_proj = Conf->n_rows; } if (Ctrl->D.active) { /* Must project the input image instead */ GMT_Report (API, GMT_MSG_INFORMATION, "Project the input image\n"); @@ -1059,8 +1496,9 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { gmt_free_header (API->GMT, &tmp_header); } - /* From here, use Grid_proj or Img_proj plus optional Intens_proj in making the (now) Cartesian rectangular image */ + /* From here, use Grid_proj or Img_proj plus optionally Intens_proj in making the (now) Cartesian rectangular image */ + has_content = (got_z_grid) ? false : true; /* Images always have content but grids may be all NaN */ if (use_intensity_grid) IH = gmt_get_H_hidden (Intens_proj->header); if (got_z_grid) { /* Dealing with a projected grid, so we only have one band [z]*/ Grid_proj->header->n_bands = 1; @@ -1069,8 +1507,15 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { else /* Use a different reference header for the image */ header_work = Img_proj->header; /* Later when need to refer to the header, use this copy */ - n_columns = header_work->n_columns; /* To simplify code below */ - n_rows = header_work->n_rows; + /* Assign the Conf structure members */ + Conf->Grid = Grid_proj; + Conf->Image = Img_proj; + Conf->Intens = Intens_proj; + Conf->H = header_work; + Conf->n_columns = header_work->n_columns; + Conf->n_rows = header_work->n_rows; + Conf->int_mode = use_intensity_grid; + Conf->nm = header_work->nm; /* Get or calculate a color palette file */ @@ -1083,11 +1528,15 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { } if (cpt) gmt_M_str_free (cpt); gray_only = (P && P->is_gray); /* Flag that we are doing a gray scale image below */ + Conf->P = P; + if (P && P->has_pattern) GMT_Report (API, GMT_MSG_WARNING, "Patterns in CPTs will be ignored\n"); + } + if (Ctrl->W.active) { /* Check if there are just NaNs in the grid */ + for (node = 0; !has_content && node < header_work->size; node++) + if (!gmt_M_is_dnan (Conf->Grid->data[node])) has_content = true; } } - if (P && P->has_pattern) GMT_Report (API, GMT_MSG_WARNING, "Patterns in CPTs will be ignored\n"); - NaN_rgb = (P) ? P->bfn[GMT_NAN].rgb : GMT->current.setting.color_patch[GMT_NAN]; /* Determine which color represents a NaN grid node */ if (Ctrl->Q.active) { /* Want colormasking via the grid's NaN entries */ if (gray_only) { @@ -1113,8 +1562,8 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { img_wesn[YHI] = GMT->current.proj.rect_m[YHI]; img_wesn[YLO] = GMT->current.proj.rect_m[YLO]; } /* Determine raster image pixel sizes in the final units */ - img_inc[0] = (img_wesn[XHI] - img_wesn[XLO]) / (n_columns - !grid_registration); - img_inc[1] = (img_wesn[YHI] - img_wesn[YLO]) / (n_rows - !grid_registration); + img_inc[0] = (img_wesn[XHI] - img_wesn[XLO]) / (Conf->n_columns - !grid_registration); + img_inc[1] = (img_wesn[YHI] - img_wesn[YLO]) / (Conf->n_rows - !grid_registration); if (grid_registration == GMT_GRID_NODE_REG) { /* Adjust domain by 1/2 pixel since SW pixel is outside the domain */ img_wesn[XLO] -= 0.5 * img_inc[0]; img_wesn[XHI] += 0.5 * img_inc[0]; @@ -1166,8 +1615,8 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { if (Ctrl->M.active || gray_only) /* Only need a byte-array to hold this image */ bitimage_8 = gmt_M_memory (GMT, NULL, header_work->nm, unsigned char); else { /* Need 3-byte array for a 24-bit image, plus possibly 3 bytes for the NaN mask color */ - if (Ctrl->Q.active) colormask_offset = 3; - bitimage_24 = gmt_M_memory (GMT, NULL, 3 * header_work->nm + colormask_offset, unsigned char); + if (Ctrl->Q.active) Conf->colormask_offset = 3; + bitimage_24 = gmt_M_memory (GMT, NULL, 3 * header_work->nm + Conf->colormask_offset, unsigned char); } if (P && Ctrl->Q.active && !(Ctrl->M.active || gray_only)) { for (k = 0; k < 3; k++) bitimage_24[k] = gmt_M_u255 (P->bfn[GMT_NAN].rgb[k]); /* Scale the NaN rgb up to 0-255 range */ @@ -1180,107 +1629,58 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { normal_x = !(GMT->current.proj.projection_GMT == GMT_LINEAR && !GMT->current.proj.xyz_pos[0] && !resampled); normal_y = !(GMT->current.proj.projection_GMT == GMT_LINEAR && !GMT->current.proj.xyz_pos[1] && !resampled); - actual_row = gmt_M_memory (GMT, NULL, n_rows, unsigned int); /* Deal with any reversal of the y-axis due to -J */ - for (row = 0; row < n_rows; row++) actual_row[row] = (normal_y) ? row : n_rows - row - 1; - actual_col = gmt_M_memory (GMT, NULL, n_columns, unsigned int); /* Deal with any reversal of the x-axis due to -J */ - for (col = 0; col < n_columns; col++) actual_col[col] = (normal_x) ? col : n_columns - col - 1; + Conf->actual_row = gmt_M_memory (GMT, NULL, Conf->n_rows, unsigned int); /* Deal with any reversal of the y-axis due to -J */ + for (row = 0; row < Conf->n_rows; row++) Conf->actual_row[row] = (normal_y) ? row : Conf->n_rows - row - 1; + Conf->actual_col = gmt_M_memory (GMT, NULL, Conf->n_columns, unsigned int); /* Deal with any reversal of the x-axis due to -J */ + for (col = 0; col < Conf->n_columns; col++) Conf->actual_col[col] = (normal_x) ? col : Conf->n_columns - col - 1; - intensity_mode = use_intensity_grid; /* Set bit 1 */ - if (!got_z_grid || (IH && IH->reset_pad)) intensity_mode |= 2; /* Add bit 2 */ rgb_cube_scan = (P && Ctrl->Q.active && !Ctrl->A.active); /* Need to look for unique rgb for PostScript masking */ - step = (gray_only || Ctrl->M.active) ? 1 : 3; /* Evaluate colors at least once (try = 0), but may do twice if -Q is active and we need to select another unique NaN color not used in the image */ for (try = 0, done = false; !done && try < 2; try++) { - if (got_z_grid && (!Ctrl->Q.active || Ctrl->A.active)) { /* Got a grid and need to look up color via the CPT */ - int srow, scol; /* Due to OPENMP on Windows requiring signed int loop variables */ - GMT_Report (API, GMT_MSG_INFORMATION, "Basic z(x,y) with optional illumination and no PostScript colormasking.\n"); -#ifdef _OPENMP -#pragma omp parallel for private(srow,byte,kk,scol,node,index,rgb,k) shared(n_rows,header_work,actual_row,n_columns,step,GMT,P,Grid_proj,Ctrl,intensity_mode,Intens_proj,gray_only,bitimage_8,bitimage_24) -#endif - for (srow = 0; srow < (int)n_rows; srow++) { /* March along scanlines */ - byte = colormask_offset + srow * n_columns * step; - kk = gmt_M_ijpgi (header_work, actual_row[srow], 0); /* Start pixel of this row */ - for (scol = 0; scol < (int)n_columns; scol++) { /* Compute rgb for each pixel along this scanline */ - node = kk + actual_col[scol]; - index = gmt_get_rgb_from_z (GMT, P, Grid_proj->data[node], rgb); - if (index != (GMT_NAN - 3)) has_content = true; - if (index != (GMT_NAN - 3) && Ctrl->I.active) { /* Need to deal with illumination */ - if (intensity_mode & 1) /* Intensity value comes from the grid */ - gmt_illuminate (GMT, Intens_proj->data[node], rgb); - else /* A constant (ambient) intensity was given via -I */ - gmt_illuminate (GMT, Ctrl->I.value, rgb); - } - /* Assigning bytes to colorimage */ - if (gray_only) /* Color table only has grays, just use r since r = g = b here */ - bitimage_8[byte] = gmt_M_u255 (rgb[0]); - else if (Ctrl->M.active) /* Convert rgb to gray using the gmt_M_yiq transformation */ - bitimage_8[byte] = gmt_M_u255 (gmt_M_yiq (rgb)); - else { /* Here we do r/g/b 24-bit color */ - for (k = 0; k < 3; k++) bitimage_24[byte+k] = gmt_M_u255 (rgb[k]); - } - byte += step; - } + if (got_z_grid) { /* Dealing with Grids and CPT lookup */ + if (Ctrl->I.active) { /* With intensity */ + if (gray_only) /* Grid, grayscale, with intensity, with intensity */ + grdimage_grd_gray_with_intensity (GMT, Ctrl, Conf, bitimage_8); + else if (Ctrl->M.active) /* Grid, color converted to gray, with intensity */ + grdimage_grd_c2s_with_intensity (GMT, Ctrl, Conf, bitimage_8); + else if (Ctrl->Q.active) /* Must deal with colormasking */ + grdimage_grd_color_with_intensity_CM (GMT, Ctrl, Conf, bitimage_24, rgb_used); + else /* Grid, color, with intensity, no colormasking */ + grdimage_grd_color_with_intensity (GMT, Ctrl, Conf, bitimage_24); + } + else { /* No intensity */ + if (gray_only) /* Grid, grayscale, no intensity */ + grdimage_grd_gray_no_intensity (GMT, Ctrl, Conf, bitimage_8); + else if (Ctrl->M.active) /* Grid, color converted to gray, with intensity */ + grdimage_grd_c2s_no_intensity (GMT, Ctrl, Conf, bitimage_8); + else if (Ctrl->Q.active) /* Must deal with colormasking */ + grdimage_grd_color_no_intensity_CM (GMT, Ctrl, Conf, bitimage_24, rgb_used); + else /* Grid, color, no intensity */ + grdimage_grd_color_no_intensity (GMT, Ctrl, Conf, bitimage_24); } - done = true; /* Only doing the loop once here since no -Q */ } - else { /* Dealing with images and/or PostScript colormasking */ - /* Because we fill bitimage in via scanlines, there are complications when the image has reversed x and/or y direction (normal_x, normal_y are false) */ - for (row = 0, byte = colormask_offset; row < n_rows; row++) { /* March along scanlines in the output bitimage */ - kk = gmt_M_ijpgi (header_work, actual_row[row], 0); /* Start pixel of this row */ - if (Ctrl->D.active) { /* Must update pixel node node_RGBA once or for each row */ - if ((row == 0 || !normal_y || !normal_x)) node_RGBA = kk; /* First time per row equals 'node', afterwards it grows alone, except when one of normal_? is false */ - if (!normal_x) node_RGBA += (n_columns-1)*Img_proj->header->n_bands; /* Must start at last column instead of first column */ - } - for (col = 0; col < n_columns; col++) { /* Compute rgb for each pixel along this scanline */ - node = kk + actual_col[col]; /* Index into the current grid node value */ - if (got_z_grid) { /* Got a grid and need to look up color via the CPT */ - index = gmt_get_rgb_from_z (GMT, P, Grid_proj->data[node], rgb); - if (index != (GMT_NAN - 3)) has_content = true; - } - else { /* Input was an image */ - if (gray_only) /* Just pick the next byte as gray shade */ - rgb[0] = gmt_M_is255 (Img_proj->data[node_RGBA++]); - else { /* Get the rgb triplet from the image */ - for (k = 0; k < 3; k++) rgb[k] = gmt_M_is255 (Img_proj->data[node_RGBA++]); - if (Img_proj->header->n_bands == 4) { /* Dealing with an image with transparency values */ - /* JL: Here we assume background color is white, hence t * 1. - But what would it take to have a user selected background color? */ - double o, t; /* o - opacity, t = transparency */ - o = Img_proj->data[node_RGBA] / 255.0; t = 1 - o; - rgb[0] = o * rgb[0] + t * GMT->current.map.frame.fill[GMT_Z].rgb[0]; - rgb[1] = o * rgb[1] + t * GMT->current.map.frame.fill[GMT_Z].rgb[1]; - rgb[2] = o * rgb[2] + t * GMT->current.map.frame.fill[GMT_Z].rgb[2]; - node_RGBA++; /* Skip past that transparency byte */ - } - } - if (!normal_x) node_RGBA -= 2*Img_proj->header->n_bands; /* Must go to the start of previous column instead of where we are (start of next) */ - } - - if (Ctrl->I.active && index != (GMT_NAN - 3)) { /* Need to deal with optional illumination for non-NAN pixels */ - if (intensity_mode & 1) { /* Intensity value comes from the grid */ - if (intensity_mode & 2) /* Must recompute "node" with the gmt_M_ijp macro suitable for the image */ - node = gmt_M_ijp (Intens_proj->header, actual_row[row], actual_col[col]); - gmt_illuminate (GMT, Intens_proj->data[node], rgb); /* Apply illumination to this color */ - } - else /* A constant (ambient) intensity was given via -I */ - gmt_illuminate (GMT, Ctrl->I.value, rgb); /* Apply constant illumination to this color */ - } - - if (gray_only) /* Color table only has grays, just use r since r = g = b here */ - bitimage_8[byte++] = gmt_M_u255 (rgb[0]); - else if (Ctrl->M.active) /* Convert rgb to gray using the gmt_M_yiq transformation */ - bitimage_8[byte++] = gmt_M_u255 (gmt_M_yiq (rgb)); - else { /* Here we do the full r/g/b 24-bit color */ - for (k = 0; k < 3; k++) bitimage_24[byte++] = i_rgb[k] = gmt_M_u255 (rgb[k]); /* Scale up to integer 0-255 range */ - if (rgb_cube_scan && index != (GMT_NAN - 3)) /* Keep track of all r/g/b combinations used except for NaN */ - rgb_used[(i_rgb[0]*256 + i_rgb[1])*256+i_rgb[2]] = true; - } - } - if (!got_z_grid && normal_x) node_RGBA += header_work->n_bands * (header_work->pad[XLO] + header_work->pad[XHI]); /* Increment the node index for the image row, unless reverse x dir */ + else { /* Dealing with an image, so no CPT lookup nor colormasking */ + if (Ctrl->I.active) { /* With intensity */ + if (gray_only) /* Image, grayscale, with intensity, with intensity */ + grdimage_img_gray_with_intensity (GMT, Ctrl, Conf, bitimage_8); + else if (Ctrl->M.active) /* Image, color converted to gray, with intensity */ + grdimage_img_c2s_with_intensity (GMT, Ctrl, Conf, bitimage_8); + else /* Grid, color, with intensity */ + grdimage_img_color_with_intensity (GMT, Ctrl, Conf, bitimage_24); + } + else { /* No intensity */ + if (gray_only) /* Image, grayscale, no intensity */ + grdimage_img_gray_no_intensity (GMT, Ctrl, Conf, bitimage_8); + else if (Ctrl->M.active) /* Image, color converted to gray, with intensity */ + grdimage_img_c2s_no_intensity (GMT, Ctrl, Conf, bitimage_8); + else /* Image, color, no intensity */ + grdimage_img_color_no_intensity (GMT, Ctrl, Conf, bitimage_24); } } - if (rgb_cube_scan) { /* Check that we found an unused r/g/b value so that colormasking will work as advertised */ + if (Ctrl->Q.active) { /* Fill in the RGB cube use */ + int index = 0, ks; + /* Check that we found an unused r/g/b value so that colormasking will work as advertised */ index = (gmt_M_u255(P->bfn[GMT_NAN].rgb[0])*256 + gmt_M_u255(P->bfn[GMT_NAN].rgb[1]))*256 + gmt_M_u255(P->bfn[GMT_NAN].rgb[2]); /* The index into the cube for the selected NaN color */ if (rgb_used[index]) { /* This r/g/b already appears in the image as a non-NaN color; we must find a replacement NaN color */ for (index = 0, ks = -1; ks == -1 && index < 256*256*256; index++) /* Examine the entire cube */ @@ -1298,13 +1698,16 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { for (k = 0; k < 3; k++) P->bfn[GMT_NAN].rgb[k] = gmt_M_is255 (bitimage_24[k]); /* Set new NaN color */ } } + else + done = true; /* Original NaN color was never used by other nodes */ } - else /* No colormasking requested so we are done after the first pass */ - done = true; + else + done = true; /* Only doing the loop once here since no -Q */ } + if (Ctrl->Q.active) gmt_M_free (GMT, rgb_used); /* Done using the r/g/b cube */ - gmt_M_free (GMT, actual_row); - gmt_M_free (GMT, actual_col); + gmt_M_free (GMT, Conf->actual_row); + gmt_M_free (GMT, Conf->actual_col); if (use_intensity_grid) { /* Also done with the intensity grid */ if (need_to_project || !got_z_grid) { @@ -1314,8 +1717,8 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { } /* Get actual plot size of each pixel */ - dx = gmt_M_get_inc (GMT, header_work->wesn[XLO], header_work->wesn[XHI], n_columns, header_work->registration); - dy = gmt_M_get_inc (GMT, header_work->wesn[YLO], header_work->wesn[YHI], n_rows, header_work->registration); + dx = gmt_M_get_inc (GMT, header_work->wesn[XLO], header_work->wesn[XHI], Conf->n_columns, header_work->registration); + dy = gmt_M_get_inc (GMT, header_work->wesn[YLO], header_work->wesn[YHI], Conf->n_rows, header_work->registration); /* Set lower left position of image on map; this depends on grid or image registration */ @@ -1326,57 +1729,16 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { } /* Full rectangular dimension of the projected image in inches */ - x_side = dx * n_columns; - y_side = dy * n_rows; + x_side = dx * Conf->n_columns; + y_side = dy * Conf->n_rows; if (P && gray_only && !Ctrl->A.active) /* Determine if the gray image in fact is just black & white since the PostScript can then simplify */ for (kk = 0, P->is_bw = true; P->is_bw && kk < header_work->nm; kk++) if (!(bitimage_8[kk] == 0 || bitimage_8[kk] == 255)) P->is_bw = false; - if (P && P->is_bw && !Ctrl->A.active) { /* Can get away with a 1-bit image, but we must pack the original byte to 8 image bits, unless we are returning the image (8 or 24 bit only allowed) */ - int nx8, shift, b_or_w, nx_pixels; - uint64_t imsize, k8; - unsigned char *bit = NULL; - - GMT_Report (API, GMT_MSG_INFORMATION, "Converting 8-bit image to 1-bit B/W image then plotting it\n"); - - nx8 = irint (ceil (n_columns / 8.0)); /* Image width must be a multiple of 8 bits, so we round up */ - nx_pixels = nx8 * 8; /* The row length in bits after the rounding up */ - imsize = gmt_M_get_nm (GMT, nx8, n_rows); - bit = gmt_M_memory (GMT, NULL, imsize, unsigned char); /* Memory to hold the 1-bit image */ - - /* Reprocess the byte image. Here there are no worries about direction of rows, cols since that was dealt with during color assignment */ - - for (row = 0, k = k8 = 0; row < n_rows; row++) { /* Process each scan line */ - shift = 0; byte = 0; - for (col = 0; col < n_columns; col++, k++) { /* Visit each byte in the original gray shade image */ - b_or_w = (bitimage_8[k] == 255); /* Let white == 1, black == 0 */ - byte |= b_or_w; /* Add in the current bit (0 or 1) */ - shift++; /* Position us for the next bit */ - if (shift == 8) { /* Did all 8, time to dump out another byte ["another one bytes the dust", he, he,...] */ - bit[k8++] = (unsigned char) byte; /* Place the last 8 bits in output array... */ - byte = shift = 0; /* ...and start over for next sequence of 8 bit nodes */ - } - else /* Move the bits we have so far 1 step to the left */ - byte <<= 1; - } - if (shift) { /* Set the remaining bits in this bit to white; this can apply to the last 1-7 nodes on each row */ - byte |= 1; - shift++; - while (shift < 8) { - byte <<= 1; - byte |= 1; - shift++; - } - bit[k8++] = (unsigned char) byte; /* Copy final byte from this row into the image */ - } - } - - x_side = nx_pixels * dx; /* Since the image may be 1-7 bits wider than what we need we must enlarge it by a tiny amount */ - PSL_plotbitimage (PSL, x0, y0, x_side, y_side, PSL_BL, bit, nx_pixels, n_rows, Ctrl->G.rgb[GMT_FGD], Ctrl->G.rgb[GMT_BGD]); - gmt_M_free (GMT, bit); /* Done with the B/W image array */ - } - else if (gray_only || Ctrl->M.active) { /* Here we have a 1-layer 8 bit grayscale image */ + if (P && P->is_bw && !Ctrl->A.active) /* Can get away with a 1-bit image, but we must pack the original byte to 8 image bits, unless we are returning the image (8 or 24 bit only allowed) */ + grdimage_blackwhite_PS_image (GMT, Ctrl, bitimage_8, Conf->n_columns, Conf->n_rows, x0, y0, dx, dy); + else if (gray_only || Ctrl->M.active) { /* Here we have a 1-layer 8 bit grayshade image */ if (Ctrl->A.active) { /* Creating a raster image, not PostScript */ GMT_Report (API, GMT_MSG_INFORMATION, "Writing 8-bit grayshade image %s\n", way[Ctrl->A.way]); if (GMT_Write_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Out.file, Out) != GMT_NOERROR) @@ -1384,16 +1746,16 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { } else { /* Lay down a PostScript 8-bit image */ GMT_Report (API, GMT_MSG_INFORMATION, "Plotting 8-bit grayshade image\n"); - PSL_plotcolorimage (PSL, x0, y0, x_side, y_side, PSL_BL, bitimage_8, n_columns, n_rows, (Ctrl->E.device_dpi ? -8 : 8)); + PSL_plotcolorimage (PSL, x0, y0, x_side, y_side, PSL_BL, bitimage_8, Conf->n_columns, Conf->n_rows, (Ctrl->E.device_dpi ? -8 : 8)); } } else { /* Dealing with a 24-bit color image */ if (Ctrl->A.active) { /* Creating a raster image, not PostScript */ if (Ctrl->Q.active) { /* Must initialize the transparency byte (alpha): 255 everywhere except at NaNs where it should be 0 */ memset (Out->alpha, 255, header_work->nm); - for (node = row = 0; row < n_rows; row++) { + for (node = row = 0; row < Conf->n_rows; row++) { kk = gmt_M_ijpgi (header_work, row, 0); - for (col = 0; col < n_columns; col++, node++) { + for (col = 0; col < Conf->n_columns; col++, node++) { if (gmt_M_is_fnan (Grid_proj->data[kk + col])) Out->alpha[node] = 0; } } @@ -1405,10 +1767,12 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { else { /* Lay down a PostScript 24-bit color image */ GMT_Report (API, GMT_MSG_INFORMATION, "Plotting 24-bit color image\n"); PSL_plotcolorimage (PSL, x0, y0, x_side, y_side, PSL_BL, bitimage_24, (Ctrl->Q.active ? -1 : 1) * - n_columns, n_rows, (Ctrl->E.device_dpi ? -24 : 24)); + Conf->n_columns, Conf->n_rows, (Ctrl->E.device_dpi ? -24 : 24)); } } + gmt_M_free (GMT, Conf); /* Done with the configuration structure */ + if (!Ctrl->A.active) { /* Finalize PostScript plot, possibly including basemap */ if (!Ctrl->N.active) gmt_map_clip_off (GMT); @@ -1456,7 +1820,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) { Return (API->error); } - /* May return a flag that the image/PS has no data (see -Qn), or just NO_ERROR */ + /* May return a flag that the image/PS had no data (see -W), or just NO_ERROR */ ret_val = (Ctrl->W.active && !has_content) ? GMT_IMAGE_NO_DATA : GMT_NOERROR; Return (ret_val); }