/* GNUPLOT - plot3d.c */

/*[

* Copyright 1986 - 1993, 1998, 2004 Thomas Williams, Colin Kelley

*

* Permission to use, copy, and distribute this software and its

* documentation for any purpose with or without fee is hereby granted,

* provided that the above copyright notice appear in all copies and

* that both that copyright notice and this permission notice appear

* in supporting documentation.

*

* Permission to modify the software is granted, but not the right to

* distribute the complete modified source code. Modifications are to

* be distributed as patches to the released version. Permission to

* distribute binaries produced by compiling modified sources is granted,

* provided you

* 1. distribute the corresponding source modifications from the

* released version in the form of a patch file along with the binaries,

* 2. add special version identification to distinguish your version

* in addition to the base release version number,

* 3. provide your name and address as the primary contact for the

* support of your modified version, and

* 4. retain our contact information in regard to use of the base

* software.

* Permission to distribute the released version of the source code along

* with corresponding source modifications in the form of a patch file is

* granted with same provisions 2 through 4 for binary distributions.

*

* This software is provided "as is" without express or implied warranty

* to the extent permitted by applicable law.

]*/

#include "plot3d.h"

#include "gp_types.h"

#include "alloc.h"

#include "axis.h"

#include "command.h"

#include "contour.h"

#include "datafile.h"

#include "datablock.h"

#include "encoding.h"

#include "eval.h"

#include "getcolor.h"

#include "graph3d.h"

#include "hidden3d.h"

#include "interpol.h"

#include "misc.h"

#include "parse.h"

#include "pm3d.h"

#include "setshow.h"

#include "term_api.h"

#include "tabulate.h"

#include "util.h"

#include "variable.h" /* For locale handling */

#include "voxelgrid.h"

#include "plot2d.h" /* Only for store_label() */

#include "matrix.h" /* Used by thin-plate-splines in dgrid3d */

#include "help.h"

/* global variables exported by this module */

t_data_mapping mapping3d = MAP3D_CARTESIAN;

int dgrid3d_row_fineness = 10;

int dgrid3d_col_fineness = 10;

int dgrid3d_norm_value = 1;

int dgrid3d_mode = DGRID3D_QNORM;

double dgrid3d_x_scale = 1.0;

double dgrid3d_y_scale = 1.0;

TBOOLEAN dgrid3d = FALSE;

TBOOLEAN dgrid3d_kdensity = FALSE;

double boxdepth = 0.0;

/* static prototypes */

static void calculate_set_of_isolines(AXIS_INDEX value_axis, TBOOLEAN cross, struct iso_curve **this_iso,

AXIS_INDEX iso_axis, double iso_min, double iso_step, int num_iso_to_use,

AXIS_INDEX sam_axis, double sam_min, double sam_step, int num_sam_to_use

);

static int get_3ddata(struct surface_points * this_plot);

static void eval_3dplots(void);

static void grid_nongrid_data(struct surface_points * this_plot);

static void parametric_3dfixup(struct surface_points * start_plot, int *plot_num);

static struct surface_points * sp_alloc(int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2);

static void sp_replace(struct surface_points *sp, int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2);

static struct iso_curve * iso_alloc(int num);

static void iso_extend(struct iso_curve *ip, int num);

static void iso_free(struct iso_curve *ip);

/* helper functions for grid_nongrid_data() */

static double splines_kernel(double h);

static void thin_plate_splines_setup( struct iso_curve *old_iso_crvs, double **p_xx, int *p_numpoints );

static double qnorm( double dist_x, double dist_y, int q );

static double pythag( double dx, double dy );

/* helper function to detect empty data sets */

static void count_3dpoints(struct surface_points *plot, int *nt, int *ni, int *nu);

/* helper functions for parsing */

static void load_contour_label_options(struct text_label *contour_label);

typedef enum splot_component {

SP_FUNCTION, SP_DATAFILE, SP_DATABLOCK, SP_KEYENTRY, SP_VOXELGRID

} splot_component;

/* the curves/surfaces of the plot */

struct surface_points *first_3dplot = NULL;

static struct udft_entry plot_func;

int plot3d_num=0;

/* FIXME:

* Because this is global, it gets clobbered if there is more than

* one unbounded iteration in the splot command, e.g.

* splot for [i=0:*] A index i, for [j=0:*] B index j

* Moving it into (struct surface_points) would be nice but would require

* extra bookkeeping to track which plot header it is stored in.

*/

static int last_iteration_in_first_pass = INT_MAX;

/*

* It is a common mistake to try to plot a complex-valued function

* without reducing it to some derived real value like abs(f(z)).

* When this happens the user gets a not entirely obvious error message

* "All point z value undefined". Try to distinguish this by counting

* complex values as we go.

*/

static int n_complex_values = 0;

/*

* sp_alloc() allocates a surface_points structure that can hold 'num_iso_1'

* iso-curves with 'num_samp_2' samples and 'num_iso_2' iso-curves with

* 'num_samp_1' samples.

* If, however num_iso_2 or num_samp_1 is zero no iso curves are allocated.

*/

static struct surface_points *

sp_alloc(int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2)

{

struct lp_style_type default_lp_properties = DEFAULT_LP_STYLE_TYPE;

struct surface_points *sp = gp_alloc(sizeof(*sp), "surface");

memset(sp,0,sizeof(struct surface_points));

/* Initialize various fields */

sp->lp_properties = default_lp_properties;

sp->fill_properties = default_fillstyle;

if (sp->fill_properties.fillstyle == FS_EMPTY)

sp->fill_properties.fillstyle = FS_SOLID;

default_arrow_style(&(sp->arrow_properties));

if (num_iso_2 > 0 && num_samp_1 > 0) {

int i;

struct iso_curve *icrv;

for (i = 0; i < num_iso_1; i++) {

icrv = iso_alloc(num_samp_2);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

for (i = 0; i < num_iso_2; i++) {

icrv = iso_alloc(num_samp_1);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

}

return (sp);

}

/*

* sp_replace() updates a surface_points structure so it can hold 'num_iso_1'

* iso-curves with 'num_samp_2' samples and 'num_iso_2' iso-curves with

* 'num_samp_1' samples.

* If, however num_iso_2 or num_samp_1 is zero no iso curves are allocated.

*/

static void

sp_replace(

struct surface_points *sp,

int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2)

{

int i;

struct iso_curve *icrv, *icrvs = sp->iso_crvs;

while (icrvs) {

icrv = icrvs;

icrvs = icrvs->next;

iso_free(icrv);

}

sp->iso_crvs = NULL;

if (num_iso_2 > 0 && num_samp_1 > 0) {

for (i = 0; i < num_iso_1; i++) {

icrv = iso_alloc(num_samp_2);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

for (i = 0; i < num_iso_2; i++) {

icrv = iso_alloc(num_samp_1);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

} else

sp->iso_crvs = NULL;

}

/*

* sp_free() releases any memory which was previously malloc()'d to hold

* surface points.

*/

/* HBB 20000506: don't risk stack havoc by recursion, use iterative list

* cleanup instead */

void

sp_free(struct surface_points *sp)

{

while (sp) {

struct surface_points *next = sp->next_sp;

free(sp->title);

free(sp->title_position);

sp->title_position = NULL;

while (sp->contours) {

struct gnuplot_contours *next_cntrs = sp->contours->next;

free(sp->contours->coords);

free(sp->contours);

sp->contours = next_cntrs;

}

while (sp->iso_crvs) {

struct iso_curve *next_icrvs = sp->iso_crvs->next;

iso_free(sp->iso_crvs);

sp->iso_crvs = next_icrvs;

}

if (sp->labels) {

free_labels(sp->labels);

sp->labels = NULL;

}

free(sp);

sp = next;

}

}

/*

* iso_alloc() allocates a iso_curve structure that can hold 'num'

* points.

*/

struct iso_curve *

iso_alloc(int num)

{

struct iso_curve *ip;

ip = gp_alloc(sizeof(struct iso_curve), "iso curve");

ip->p_max = (num >= 0 ? num : 0);

ip->p_count = 0;

if (num > 0) {

ip->points = gp_alloc(num * sizeof(struct coordinate), "iso curve points");

memset(ip->points, 0, num * sizeof(struct coordinate));

} else

ip->points = (struct coordinate *) NULL;

ip->next = NULL;

return (ip);

}

/*

* iso_extend() reallocates a iso_curve structure to hold "num"

* points. This will either expand or shrink the storage.

*/

void

iso_extend(struct iso_curve *ip, int num)

{

if (num == ip->p_max)

return;

if (num > 0) {

ip->points = gp_realloc(ip->points, num * sizeof(struct coordinate),

"expanding 3D points");

if (num > ip->p_max)

memset( &(ip->points[ip->p_max]), 0, (num - ip->p_max) * sizeof(struct coordinate));

ip->p_max = num;

} else {

free(ip->points);

ip->points = NULL;

ip->p_max = 0;

}

}

/*

* iso_free() releases any memory which was previously malloc()'d to hold

* iso curve points.

*/

void

iso_free(struct iso_curve *ip)

{

if (ip)

free(ip->points);

free(ip);

}

void

plot3drequest()

/*

* in the parametric case we would say splot [u= -Pi:Pi] [v= 0:2*Pi] [-1:1]

* [-1:1] [-1:1] sin(v)*cos(u),sin(v)*cos(u),sin(u) in the non-parametric

* case we would say only splot [x= -2:2] [y= -5:5] sin(x)*cos(y)

*

*/

{

int dummy_token0 = -1, dummy_token1 = -1;

AXIS_INDEX axis, u_axis, v_axis;

is_3d_plot = TRUE;

if (parametric && strcmp(set_dummy_var[0], "t") == 0) {

strcpy(set_dummy_var[0], "u");

strcpy(set_dummy_var[1], "v");

}

/* initialize the arrays from the 'set' scalars */

axis_init(&axis_array[FIRST_X_AXIS], FALSE);

axis_init(&axis_array[FIRST_Y_AXIS], FALSE);

axis_init(&axis_array[FIRST_Z_AXIS], TRUE);

axis_init(&axis_array[U_AXIS], FALSE);

axis_init(&axis_array[V_AXIS], FALSE);

axis_init(&axis_array[COLOR_AXIS], TRUE);

/* Always be prepared to restore the autoscaled values on "refresh"

* Dima Kogan April 2018

* FIXME: Could we fold this into axis_init?

*/

for (axis = 0; axis < NUMBER_OF_MAIN_VISIBLE_AXES; axis++) {

AXIS *this_axis = &axis_array[axis];

if (this_axis->set_autoscale != AUTOSCALE_NONE)

this_axis->range_flags |= RANGE_WRITEBACK;

}

/* Nonlinear mapping of x or y via linkage to a hidden primary axis. */

/* The user set autoscale for the visible axis; apply it also to the hidden axis. */

for (axis = 0; axis < NUMBER_OF_MAIN_VISIBLE_AXES; axis++) {

AXIS *secondary = &axis_array[axis];

if (axis == SAMPLE_AXIS)

continue;

if (secondary->linked_to_primary

&& secondary->linked_to_primary->index == -secondary->index) {

AXIS *primary = secondary->linked_to_primary;

primary->set_autoscale = secondary->set_autoscale;

axis_init(primary, 1);

}

}

if (!term) /* unknown */

int_error(c_token, "use 'set term' to set terminal type first");

/* Range limits for the entire plot are optional but must be given */

/* in a fixed order. The keyword 'sample' terminates range parsing. */

u_axis = (parametric ? U_AXIS : FIRST_X_AXIS);

v_axis = (parametric ? V_AXIS : FIRST_Y_AXIS);

dummy_token0 = parse_range(u_axis);

dummy_token1 = parse_range(v_axis);

if (parametric) {

parse_range(FIRST_X_AXIS);

parse_range(FIRST_Y_AXIS);

}

parse_range(FIRST_Z_AXIS);

check_axis_reversed(FIRST_X_AXIS);

check_axis_reversed(FIRST_Y_AXIS);

check_axis_reversed(FIRST_Z_AXIS);

if (equals(c_token,"sample") && equals(c_token+1,"["))

c_token++;

/* Clear out any tick labels read from data files in previous plot */

for (axis=0; axis<AXIS_ARRAY_SIZE; axis++) {

struct ticdef *ticdef = &axis_array[axis].ticdef;

if (ticdef->def.user)

ticdef->def.user = prune_dataticks(ticdef->def.user);

if (!ticdef->def.user && ticdef->type == TIC_USER)

ticdef->type = TIC_COMPUTED;

}

/* use the default dummy variable unless changed */

if (dummy_token0 > 0)

copy_str(c_dummy_var[0], dummy_token0, MAX_ID_LEN);

else

strcpy(c_dummy_var[0], set_dummy_var[0]);

if (dummy_token1 > 0)

copy_str(c_dummy_var[1], dummy_token1, MAX_ID_LEN);

else

strcpy(c_dummy_var[1], set_dummy_var[1]);

/* In "set view map" mode the x2 and y2 axes are legal */

/* but must be linked to the respective primary axis. */

if (splot_map) {

if ((axis_array[SECOND_X_AXIS].ticmode && !axis_array[SECOND_X_AXIS].linked_to_primary)

|| (axis_array[SECOND_Y_AXIS].ticmode && !axis_array[SECOND_Y_AXIS].linked_to_primary))

int_error(NO_CARET,

"Secondary axis must be linked to primary axis in order to draw tics");

}

eval_3dplots();

}

/* Helper function for refresh command. Reexamine each data point and update the

* flags for INRANGE/OUTRANGE/UNDEFINED based on the current limits for that axis.

* Normally the axis limits are already known at this point. But if the user has

* forced "set autoscale" since the previous plot or refresh, we need to reset the

* axis limits and try to approximate the full auto-scaling behaviour.

*/

void

refresh_3dbounds(struct surface_points *first_plot, int nplots)

{

struct surface_points *this_plot = first_plot;

int iplot; /* plot index */

for (iplot = 0; iplot < nplots; iplot++, this_plot = this_plot->next_sp) {

int i; /* point index */

struct axis *x_axis = &axis_array[FIRST_X_AXIS];

struct axis *y_axis = &axis_array[FIRST_Y_AXIS];

struct axis *z_axis = &axis_array[FIRST_Z_AXIS];

struct iso_curve *this_curve;

/* IMAGE clipping is done elsewhere, so we don't need INRANGE/OUTRANGE checks */

if (this_plot->plot_style == IMAGE

|| this_plot->plot_style == RGBIMAGE

|| this_plot->plot_style == RGBA_IMAGE) {

if (x_axis->set_autoscale)

process_image(this_plot, IMG_UPDATE_AXES);

continue;

}

for (this_curve = this_plot->iso_crvs; this_curve; this_curve = this_curve->next) {

/* VECTOR plots consume two iso_crvs structures, one for heads and one for tails.

* Only the first one has the true point count; the second one claims zero.

* FIXME: Maybe we should change this?

*/

int n_points;

if (this_plot->plot_style == VECTOR)

n_points = this_plot->iso_crvs->p_count;

else

n_points = this_curve->p_count;

for (i=0; i<n_points; i++) {

struct coordinate *point = &this_curve->points[i];

if (point->type == UNDEFINED)

continue;

else

point->type = INRANGE;

/* If the state has been set to autoscale since the last plot,

* mark everything INRANGE and re-evaluate the axis limits now.

* Otherwise test INRANGE/OUTRANGE against previous axis limits.

*/

/* This autoscaling logic is parallel to that in

* refresh_bounds() in plot2d.c

*/

if (!this_plot->noautoscale) {

autoscale_one_point(x_axis, point->x);

autoscale_one_point(y_axis, point->y);

}

if (!inrange(point->x, x_axis->min, x_axis->max)

|| !inrange(point->y, y_axis->min, y_axis->max)) {

point->type = OUTRANGE;

continue;

}

if (!this_plot->noautoscale) {

autoscale_one_point(z_axis, point->z);

}

if (!inrange(point->z, z_axis->min, z_axis->max)) {

point->type = OUTRANGE;

continue;

}

} /* End of points in this curve */

} /* End of curves in this plot */

} /* End of plots in this splot command */

/* handle 'reverse' ranges */

axis_check_range(FIRST_X_AXIS);

axis_check_range(FIRST_Y_AXIS);

axis_check_range(FIRST_Z_AXIS);

/* Make sure the bounds are reasonable, and tweak them if they aren't */

axis_checked_extend_empty_range(FIRST_X_AXIS, NULL);

axis_checked_extend_empty_range(FIRST_Y_AXIS, NULL);

axis_checked_extend_empty_range(FIRST_Z_AXIS, NULL);

}

static double

splines_kernel(double h)

{

if (h > 0.0) { return h * h * log(h); }

return 0.0;

}

/* PKJ:

This function has been hived off out of the original grid_nongrid_data().

No changes have been made, but variables only needed locally have moved

out of grid_nongrid_data() into this function. */

static void

thin_plate_splines_setup( struct iso_curve *old_iso_crvs,

double **p_xx, int *p_numpoints )

{

int i, j, k;

double *xx, *yy, *zz, *b, **K, d;

int numpoints, *indx;

struct iso_curve *oicrv;

numpoints = 0;

for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {

numpoints += oicrv->p_count;

}

xx = gp_alloc(sizeof(xx[0]) * (numpoints + 3) * (numpoints + 8),

"thin plate splines in dgrid3d");

/* the memory needed is not really (n+3)*(n+8) for now,

but might be if I take into account errors ... */

K = gp_alloc(sizeof(K[0]) * (numpoints + 3),

"matrix : thin plate splines 2d");

yy = xx + numpoints;

zz = yy + numpoints;

b = zz + numpoints;

/* HBB 20010424: Count actual input points without the UNDEFINED

* ones, as we copy them */

numpoints = 0;

for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {

struct coordinate *opoints = oicrv->points;

for (k = 0; k < oicrv->p_count; k++, opoints++) {

/* HBB 20010424: avoid crashing for undefined input */

if (opoints->type == UNDEFINED)

continue;

xx[numpoints] = opoints->x;

yy[numpoints] = opoints->y;

zz[numpoints] = opoints->z;

numpoints++;

}

}

for (i = 0; i < numpoints + 3; i++) {

K[i] = b + (numpoints + 3) * (i + 1);

}

for (i = 0; i < numpoints; i++) {

for (j = i + 1; j < numpoints; j++) {

double dx = xx[i] - xx[j], dy = yy[i] - yy[j];

K[i][j] = K[j][i] = -splines_kernel(sqrt(dx * dx + dy * dy));

}

K[i][i] = 0.0; /* here will come the weights for errors */

b[i] = zz[i];

}

for (i = 0; i < numpoints; i++) {

K[i][numpoints] = K[numpoints][i] = 1.0;

K[i][numpoints + 1] = K[numpoints + 1][i] = xx[i];

K[i][numpoints + 2] = K[numpoints + 2][i] = yy[i];

}

b[numpoints] = 0.0;

b[numpoints + 1] = 0.0;

b[numpoints + 2] = 0.0;

K[numpoints][numpoints] = 0.0;

K[numpoints][numpoints + 1] = 0.0;

K[numpoints][numpoints + 2] = 0.0;

K[numpoints + 1][numpoints] = 0.0;

K[numpoints + 1][numpoints + 1] = 0.0;

K[numpoints + 1][numpoints + 2] = 0.0;

K[numpoints + 2][numpoints] = 0.0;

K[numpoints + 2][numpoints + 1] = 0.0;

K[numpoints + 2][numpoints + 2] = 0.0;

indx = gp_alloc(sizeof(indx[0]) * (numpoints + 3), "indexes lu");

/* actually, K is *not* positive definite, but

has only non zero real eigenvalues ->

we can use an lu_decomp safely */

lu_decomp(K, numpoints + 3, indx, &d);

lu_backsubst(K, numpoints + 3, indx, b);

free( K );

free( indx );

*p_xx = xx;

*p_numpoints = numpoints;

}

static double

qnorm( double dist_x, double dist_y, int q )

{

double dist = 0.0;

#if (1)

switch (q) {

case 1:

dist = pythag(dist_x, dist_y);

break;

case 2:

dist = dist_x*dist_x + dist_y*dist_y;

break;

default:

dist = pow( pythag(dist_x, dist_y), q);

break;

}

#else /* old code (versions 3.5 - 5.2) */

switch (q) {

case 1:

dist = dist_x + dist_y;

break;

case 2:

dist = dist_x * dist_x + dist_y * dist_y;

break;

case 4:

dist = dist_x * dist_x + dist_y * dist_y;

dist *= dist;

break;

case 8:

dist = dist_x * dist_x + dist_y * dist_y;

dist *= dist;

dist *= dist;

break;

case 16:

dist = dist_x * dist_x + dist_y * dist_y;

dist *= dist;

dist *= dist;

dist *= dist;

break;

default:

dist = pow(dist_x, (double)q ) + pow(dist_y, (double)q );

break;

}

#endif

return dist;

}

/* This is from Numerical Recipes in C, 2nd ed, p70 */

static double

pythag( double dx, double dy )

{

double x, y;

x = fabs(dx);

y = fabs(dy);

if (x > y) { return x*sqrt(1.0 + (y*y)/(x*x)); }

if (y==0.0) { return 0.0; }

return y*sqrt(1.0 + (x*x)/(y*y));

}

static void

grid_nongrid_data(struct surface_points *this_plot)

{

int i, j, k;

double x, y, z, w, dx, dy, xmin, xmax, ymin, ymax;

double c;

coord_type dummy_type;

struct iso_curve *old_iso_crvs = this_plot->iso_crvs;

struct iso_curve *icrv, *oicrv, *oicrvs;

/* these are only needed for thin_plate_splines */

double *xx = NULL, *yy = NULL, *zz = NULL, *b = NULL;

int numpoints = 0;

/* Version 5.3 - allow gridding of separate color column so long as

* we don't need to generate splines for it

*/

if (this_plot->pm3d_color_from_column && dgrid3d_mode == DGRID3D_SPLINES)

int_error(NO_CARET, "Spline gridding of a separate color column is not implemented");

/* Compute XY bounding box on the original data. */

/* FIXME HBB 20010424: Does this make any sense? Shouldn't we just

* use whatever the x and y ranges have been found to be, and

* that's that? The largest difference this is going to make is if

* we plot a datafile that doesn't span the whole x/y range

* used. Do we want a dgrid3d over the actual data rectangle, or

* over the xrange/yrange area? */

xmin = xmax = old_iso_crvs->points[0].x;

ymin = ymax = old_iso_crvs->points[0].y;

for (icrv = old_iso_crvs; icrv != NULL; icrv = icrv->next) {

struct coordinate *points = icrv->points;

for (i = 0; i < icrv->p_count; i++, points++) {

/* HBB 20010424: avoid crashing for undefined input */

if (points->type == UNDEFINED)

continue;

if (xmin > points->x)

xmin = points->x;

if (xmax < points->x)

xmax = points->x;

if (ymin > points->y)

ymin = points->y;

if (ymax < points->y)

ymax = points->y;

}

}

dx = (xmax - xmin) / (dgrid3d_col_fineness - 1);

dy = (ymax - ymin) / (dgrid3d_row_fineness - 1);

/* Create the new grid structure, and compute the low pass filtering from

* non grid to grid structure.

*/

this_plot->iso_crvs = NULL;

this_plot->num_iso_read = dgrid3d_col_fineness;

this_plot->has_grid_topology = TRUE;

if (dgrid3d_mode == DGRID3D_SPLINES) {

thin_plate_splines_setup( old_iso_crvs, &xx, &numpoints );

yy = xx + numpoints;

zz = yy + numpoints;

b = zz + numpoints;

}

for (i = 0, x = xmin; i < dgrid3d_col_fineness; i++, x += dx) {

struct coordinate *points;

icrv = iso_alloc(dgrid3d_row_fineness + 1);

icrv->p_count = dgrid3d_row_fineness;

icrv->next = this_plot->iso_crvs;

this_plot->iso_crvs = icrv;

points = icrv->points;

for (j=0, y=ymin; j<dgrid3d_row_fineness; j++, y+=dy, points++) {

z = w = 0.0;

c = 0.0;

/* as soon as ->type is changed to UNDEFINED, break out of

* two inner loops! */

points->type = INRANGE;

if (dgrid3d_mode == DGRID3D_SPLINES) {

z = b[numpoints];

for (k = 0; k < numpoints; k++) {

double dx = xx[k] - x, dy = yy[k] - y;

z = z - b[k] * splines_kernel(sqrt(dx * dx + dy * dy));

}

z = z + b[numpoints + 1] * x + b[numpoints + 2] * y;

} else { /* everything, except splines */

for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {

struct coordinate *opoints = oicrv->points;

for (k = 0; k < oicrv->p_count; k++, opoints++) {

if (dgrid3d_mode == DGRID3D_QNORM) {

double dist = qnorm( fabs(opoints->x - x),

fabs(opoints->y - y),

dgrid3d_norm_value );

if (dist == 0.0) {

/* HBB 981209: flag the first undefined z and return */

points->type = UNDEFINED;

z = opoints->z;

c = opoints->CRD_COLOR;

w = 1.0;

break; /* out of inner loop */

} else {

z += opoints->z / dist;

c += opoints->CRD_COLOR / dist;

w += 1.0/dist;

}

} else { /* ALL else: not spline, not qnorm! */

double weight = 0.0;

double dist=pythag((opoints->x-x)/dgrid3d_x_scale,

(opoints->y-y)/dgrid3d_y_scale);

if (dgrid3d_mode == DGRID3D_GAUSS) {

weight = exp( -dist*dist );

} else if (dgrid3d_mode == DGRID3D_CAUCHY) {

weight = 1.0/(1.0 + dist*dist );

} else if (dgrid3d_mode == DGRID3D_EXP) {

weight = exp( -dist );

} else if (dgrid3d_mode == DGRID3D_BOX) {

weight = (dist<1.0) ? 1.0 : 0.0;

} else if (dgrid3d_mode == DGRID3D_HANN) {

if (dist < 1.0)

weight = 0.5*(1+cos(M_PI*dist));

}

z += opoints->z * weight;

c += opoints->CRD_COLOR * weight;

w += weight;

}

}

/* PKJ: I think this is only relevant for qnorm */

if (points->type != INRANGE)

break; /* out of the second-inner loop as well ... */

}

} /* endif( dgrid3d_mode == DGRID3D_SPLINES ) */

/* Now that we've escaped the loops safely, we know that we

* do have a good value in z and w, so we can proceed just as

* if nothing had happened at all. Nice, isn't it? */

points->type = INRANGE;

points->x = x;

points->y = y;

/* Honor requested x and y limits */

/* Historical note: This code was not in 4.0 or 4.2. It imperfectly */

/* restores the clipping behaviour of version 3.7 and earlier. */

if ((x < axis_array[x_axis].min && !(axis_array[x_axis].autoscale & AUTOSCALE_MIN))

|| (x > axis_array[x_axis].max && !(axis_array[x_axis].autoscale & AUTOSCALE_MAX))

|| (y < axis_array[y_axis].min && !(axis_array[y_axis].autoscale & AUTOSCALE_MIN))

|| (y > axis_array[y_axis].max && !(axis_array[y_axis].autoscale & AUTOSCALE_MAX)))

points->type = OUTRANGE;

if (dgrid3d_mode != DGRID3D_SPLINES && !dgrid3d_kdensity) {

z = z / w;

c = c / w;

}

STORE_AND_UPDATE_RANGE(points->z, z, points->type, z_axis,

this_plot->noautoscale, continue);

if (!this_plot->pm3d_color_from_column)

c = z;

dummy_type = points->type;

STORE_AND_UPDATE_RANGE( points->CRD_COLOR, c, dummy_type, COLOR_AXIS,

this_plot->noautoscale, continue);

}

}

free(xx); /* safe to call free on NULL pointer if splines not used */

/* Delete the old non grid data. */

for (oicrvs = old_iso_crvs; oicrvs != NULL;) {

oicrv = oicrvs;

oicrvs = oicrvs->next;

iso_free(oicrv);

}

}

/* Get 3D data from file, and store into this_plot data

* structure. Takes care of 'set mapping' and 'set dgrid3d'.

*

* Notice: this_plot->token is end of datafile spec, before title etc

* will be moved past title etc after we return */

static int

get_3ddata(struct surface_points *this_plot)

{

int xdatum = 0;

int ydatum = 0;

int j;

int pt_in_iso_crv = 0;

struct iso_curve *this_iso;

int retval = 0;

double v[MAXDATACOLS];

/* Initialize the space that will hold input data values */

memset(v, 0, sizeof(v));

if (mapping3d == MAP3D_CARTESIAN) {

/* do this check only, if we have PM3D / PM3D-COLUMN not compiled in */

if (df_no_use_specs == 2)

int_error(this_plot->token, "Need 1 or 3 columns for cartesian data");

/* HBB NEW 20060427: if there's only one, explicit using

* column, it's z data. df_axis[] has to reflect that, so

* df_readline() will expect time/date input. */

if (df_no_use_specs == 1)

df_axis[0] = FIRST_Z_AXIS;

} else {

if (df_no_use_specs == 1)

int_error(this_plot->token, "Need 2 or 3 columns for polar data");

}

this_plot->num_iso_read = 0;

this_plot->has_grid_topology = TRUE;

this_plot->pm3d_color_from_column = FALSE;

/* we ought to keep old memory - most likely case

* is a replot, so it will probably exactly fit into

* memory already allocated ?

*/

if (this_plot->iso_crvs != NULL) {

struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;

while (icrvs) {

icrv = icrvs;

icrvs = icrvs->next;

iso_free(icrv);

}

this_plot->iso_crvs = NULL;

}

/* data file is already open */

if (df_matrix)

this_plot->has_grid_topology = TRUE;

{

/*{{{ read surface from text file */

struct iso_curve *local_this_iso = iso_alloc(samples_1);

struct coordinate *cp;

struct coordinate *cphead = NULL; /* Only for VECTOR plots */

double x, y, z;

double xlow = 0, xhigh = 0;

double xtail, ytail, ztail;

double weight = 1.0;

double zlow = VERYLARGE, zhigh = -VERYLARGE;

double color = VERYLARGE;

int pm3d_color_from_column = FALSE;

#define color_from_column(x) pm3d_color_from_column = x

if (this_plot->plot_style == LABELPOINTS)

expect_string( 4 );

if (this_plot->plot_style == VECTOR) {

local_this_iso->next = iso_alloc(samples_1);

local_this_iso->next->p_count = 0;

}

if (this_plot->plot_style == POLYGONS) {

this_plot->has_grid_topology = FALSE;

track_pm3d_quadrangles = TRUE;

}

/* If the user has set an explicit locale for numeric input, apply it */

/* here so that it affects data fields read from the input file. */

set_numeric_locale();

/* Initial state */

df_warn_on_missing_columnheader = TRUE;

while ((retval = df_readline(v,MAXDATACOLS)) != DF_EOF) {

j = retval;

if (j == 0) /* not blank line, but df_readline couldn't parse it */

int_warn(NO_CARET, "Bad data on line %d of file %s",

df_line_number, df_filename ? df_filename : "");

if (j == DF_SECOND_BLANK)

break; /* two blank lines */

if (j == DF_FIRST_BLANK) {

/* Images are in a sense similar to isocurves.

* However, the routine for images is written to

* compute the two dimensions of coordinates by

* examining the data alone. That way it can be used

* in the 2D plots, for which there is no isoline

* record. So, toss out isoline information for

* images.

*/

if ((this_plot->plot_style == IMAGE)

|| (this_plot->plot_style == RGBIMAGE)

|| (this_plot->plot_style == RGBA_IMAGE))

continue;

if (this_plot->plot_style == VECTOR)

continue;

/* one blank line */

if (pt_in_iso_crv == 0) {

if (xdatum == 0)

continue;

pt_in_iso_crv = xdatum;

}

if (xdatum > 0) {

local_this_iso->p_count = xdatum;

local_this_iso->next = this_plot->iso_crvs;

this_plot->iso_crvs = local_this_iso;

this_plot->num_iso_read++;

if (xdatum != pt_in_iso_crv)

this_plot->has_grid_topology = FALSE;

local_this_iso = iso_alloc(pt_in_iso_crv);

xdatum = 0;

ydatum++;

}

continue;

}

else if (j == DF_FOUND_KEY_TITLE){

/* only the shared part of the 2D and 3D headers is used */

df_set_key_title((struct curve_points *)this_plot);

continue;

}

else if (j == DF_KEY_TITLE_MISSING){

fprintf(stderr,

"get_data: key title not found in requested column\n"

);

continue;

}

else if (j == DF_COLUMN_HEADERS) {

continue;

}

/* its a data point or undefined */

if (xdatum >= local_this_iso->p_max) {

/* overflow about to occur. Extend size of points[]

* array. Double the size, and add 1000 points, to

* avoid needlessly small steps. */

iso_extend(local_this_iso, xdatum + xdatum + 1000);