/* GNUPLOT - filters.c */
/*[
* Copyright Ethan A Merritt 2013 - 2023
* All code in this file is dual-licensed.
*
* Gnuplot license:
*
* 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.
*
* Alternative license:
*
* As an alternative to distributing code in this file under the gnuplot license,
* you may instead comply with the terms below. In this case, redistribution and
* use in source and binary forms, with or without modification, are permitted
* provided that the following conditions are met:
*
* Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer. Redistributions in binary
* form must reproduce the above copyright notice, this list of conditions and
* the following disclaimer in the documentation and/or other materials provided
* with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*
]*/
#include "filters.h"
#include "interpol.h"
#include "alloc.h"
#include "datafile.h" /* for blank_data_line */
#include "plot2d.h" /* for cp_extend() */
#include "watch.h" /* for bisect_target() */
/*
* local prototypes
*/
static int do_curve_cleanup(struct coordinate *point, int npoints);
static void winnow_interior_points (struct curve_points *plot, t_cluster *cluster);
static double fpp_SG5(struct coordinate *p);
static void cluster_stats(struct coordinate *points, t_cluster *cluster);
/* Variables related to clustering
* These really belong somewhere else, but the only routines so
* far that use cluster properties are the hull routines,
* so keep them together for now.
* cluster_outlier_threshold sets criterion for finding outliers
* not currently settable
* chi_shape_default_fraction sets the default choice of chi_length to
* chi_shape_default_fraction * max(edgelengths)
* currently controlled by "set chi_shape fraction "
*/
double cluster_outlier_threshold = 0.0;
double chi_shape_default_fraction = 0.6;
/*
* EAM December 2013
* monotonic cubic spline using the Fritsch-Carlson algorithm
* FN Fritsch & RE Carlson (1980). "Monotone Piecewise Cubic Interpolation".
* SIAM Journal on Numerical Analysis (SIAM) 17 (2): 238246. doi:10.1137/0717021.
*/
void
mcs_interp(struct curve_points *plot)
{
/* These track the original (pre-sorted) data points */
int N = plot->p_count;
struct coordinate *p = gp_realloc(plot->points, (N+1) * sizeof(coordinate), "mcs");
int i;
/* These will track the resulting smoothed curve (>= 3X original count) */
/* Larger number of samples gives smoother curve (no surprise!) */
int Nsamp = (samples_1 > 2*N) ? samples_1 : 2*N;
int Ntot = N + Nsamp;
struct coordinate *new_points = gp_alloc((Ntot) * sizeof(coordinate), "mcs");
double xstart = GPMAX(p[0].x, X_AXIS.min);
double xend = GPMIN(p[N-1].x, X_AXIS.max);
double xstep = (xend - xstart) / (Nsamp - 1);
/* Load output x coords for sampling */
for (i=0; ip_count = Ntot;
plot->p_max = Ntot + 1;
#undef DX
#undef SLOPE
#undef C1
#undef C2
#undef C3
}
/*
* Binned histogram of input values.
*
* plot FOO using N:(1) bins{=} {binrange=[binlow:binhigh]}
* {binwidth=} with boxes
*
* If no binrange is given, binlow and binhigh are taken from the x range of the data.
* In either of these cases binlow is the midpoint x-coordinate of the first bin
* and binhigh is the midpoint x-coordinate of the last bin.
* Points that lie exactly on a bin boundary are assigned to the upper bin.
* Bin assignments are not affected by "set xrange".
* Notes:
* binwidth = (binhigh-binlow) / (nbins-1)
* xmin = binlow - binwidth/2
* xmax = binhigh + binwidth/2
* first bin holds points with (xmin =< x < xmin + binwidth)
* last bin holds points with (xmax-binwidth =< x < binhigh + binwidth)
*
* binopt = 0 (default) return sum of y values for points in each bin
* binopt = 1 return mean of y values for points in each bin
*
* Ethan A Merritt 2015
*/
void
make_bins(struct curve_points *plot, int nbins,
double binlow, double binhigh, double binwidth,
int binopt)
{
int i, binno;
double *bin;
double bottom, top, range;
int *members;
struct axis *xaxis = &axis_array[plot->x_axis];
struct axis *yaxis = &axis_array[plot->y_axis];
double ymax = 0;
int N = plot->p_count;
/* Find the range of points to be binned */
if (binlow != binhigh) {
/* Explicit binrange [min:max] in the plot command */
bottom = binlow;
top = binhigh;
} else {
/* Take binrange from the data itself */
bottom = VERYLARGE; top = -VERYLARGE;
for (i=0; i plot->points[i].x)
bottom = plot->points[i].x;
if (top < plot->points[i].x)
top = plot->points[i].x;
}
if (top 0) {
double temp;
nbins = 1 + (top - bottom) / binwidth;
temp = nbins * binwidth - (top - bottom);
bottom -= temp/2.;
top += temp/2.;
}
/* otherwise we use (N-1) intervals between midpoints of bin 1 and bin N */
else {
binwidth = (top - bottom) / (nbins - 1);
bottom -= binwidth/2.;
top += binwidth/2.;
}
range = top - bottom;
bin = gp_alloc(nbins*sizeof(double), "bins");
members = gp_alloc(nbins*sizeof(int), "bins");
for (i=0; ipoints[i].x - bottom) / range);
if (0 points[i].y;
members[binno]++;
}
}
if (xaxis->autoscale & AUTOSCALE_MIN) {
if (xaxis->min > bottom)
xaxis->min = bottom;
}
if (xaxis->autoscale & AUTOSCALE_MAX) {
if (xaxis->max < top)
xaxis->max = top;
}
/* Replace the original data with one entry per bin.
* new x = midpoint of bin
* new y = sum of individual y values over all points in bin
* new z = number of points in the bin
*/
plot->p_count = nbins;
cp_extend(plot, nbins);
for (i=0; i 1))
ybin = bin[i]/members[i];
plot->points[i].type = INRANGE;
plot->points[i].x = bincent;
plot->points[i].xlow = bincent - binwidth/2.;
plot->points[i].xhigh = bincent + binwidth/2.;
plot->points[i].y = ybin;
plot->points[i].ylow = ybin;
plot->points[i].yhigh = ybin;
plot->points[i].z = members[i];
if (inrange(bincent, xaxis->min, xaxis->max)) {
if (ymax < ybin)
ymax = ybin;
} else {
plot->points[i].type = OUTRANGE;
}
}
if (yaxis->autoscale & AUTOSCALE_MIN) {
if (yaxis->min > 0)
yaxis->min = 0;
}
if (yaxis->autoscale & AUTOSCALE_MAX) {
if (yaxis->max < ymax)
yaxis->max = ymax;
}
/* Recheck range on y */
for (i=0; ipoints[i].y, yaxis->min, yaxis->max))
plot->points[i].type = OUTRANGE;
/* Clean up */
free(bin);
free(members);
}
/*
* spline approximation of 3D lines
* do_3d_cubic gen_2d_path_splines gen_3d_splines
* Ethan A Merritt 2019
*/
/*
* Replace one isocurve with a 3D natural cubic spline interpolation.
* If there are multiple isocurves, or multiple curves with isocurves,
* the caller must sort that out and call here separately for each one.
* TODO:
* number of spline samples should be independent of "set samples"
*/
static void
do_3d_cubic(struct iso_curve *curve, enum PLOT_SMOOTH smooth_option)
{
int i, l;
int nseg = samples_1;
struct coordinate *old_points, *new_points;
double xrange, yrange, zrange;
double dx, dy, dz;
double maxdx, maxdy, maxdz;
double t, tsum, tstep;
spline_coeff *sc_x = NULL;
spline_coeff *sc_y = NULL;
spline_coeff *sc_z = NULL;
old_points = curve->points;
/*
* Sanity check axis ranges.
* This catches curves that lie in a plane of constant x or y.
* The fixup prints a warning to the user but we don't see it here.
*/
axis_checked_extend_empty_range(FIRST_X_AXIS, "at time of spline generation");
axis_checked_extend_empty_range(FIRST_Y_AXIS, "at time of spline generation");
/* prevent gross mismatch of x/y/z units */
xrange = fabs(axis_array[FIRST_X_AXIS].max - axis_array[FIRST_X_AXIS].min);
yrange = fabs(axis_array[FIRST_Y_AXIS].max - axis_array[FIRST_Y_AXIS].min);
zrange = fabs(axis_array[FIRST_Z_AXIS].max - axis_array[FIRST_Z_AXIS].min);
/* Construct path-length vector; store it in unused slot of old_points */
t = tsum = 0.0;
maxdx = maxdy = maxdz = 0.0;
old_points[0].CRD_PATH = 0;
for (i = 1; i < curve->p_count; i++) {
dx = (old_points[i].x - old_points[i-1].x) / xrange;
dy = (old_points[i].y - old_points[i-1].y) / yrange;
dz = (old_points[i].z - old_points[i-1].z) / zrange;
tsum += sqrt( dx*dx + dy*dy + dz*dz );
old_points[i].CRD_PATH = tsum;
/* Track planarity */
if (fabs(dx) > maxdx)
maxdx = fabs(dx);
if (fabs(dy) > maxdy)
maxdy = fabs(dy);
if (fabs(dz) > maxdz)
maxdz = fabs(dz);
}
/* Normalize so that the path always runs from 0 to 1 */
for (i = 1; i < curve->p_count; i++)
old_points[i].CRD_PATH /= tsum;
tstep = old_points[curve->p_count-1].CRD_PATH / (double)(nseg - 1);
/* Create new list to hold interpolated points */
new_points = gp_alloc((nseg+1) * sizeof(struct coordinate), "3D spline");
memset( new_points, 0, (nseg+1) * sizeof(struct coordinate));
/*
* If the curve being fitted lies entirely in one plane,
* we can do better by fitting a 2D spline rather than a 3D spline.
* This benefits the relatively common case of drawing a stack of
* 2D plots (e.g. fence plots).
* First check for a curve lying in the yz plane (x = constant).
*/
if (maxdx < FLT_EPSILON) {
tstep = (old_points[curve->p_count-1].y - old_points[0].y) / (double)(nseg - 1);
if (smooth_option == SMOOTH_ACSPLINES)
sc_z = cp_approx_spline(curve->points, curve->p_count, 1, 2, 3);
else
sc_z = cp_tridiag(curve->points, curve->p_count, 1, 2);
for (i = 0, l = 0; i < nseg; i++) {
double temp;
t = old_points[0].y + i * tstep;
/* Move forward to the spline interval this point is in */
while ((t >= old_points[l + 1].y) && (l < curve->p_count- 2))
l++;
temp = t - old_points[l].y;
new_points[i].x = old_points[l].x; /* All the same */
new_points[i].y = t;
new_points[i].z = ((sc_z[l][3] * temp + sc_z[l][2]) * temp + sc_z[l][1])
* temp + sc_z[l][0];
}
}
/*
* Check for a curve lying in the xz plane (y = constant).
*/
else if (maxdy < FLT_EPSILON) {
tstep = (old_points[curve->p_count-1].x - old_points[0].x) / (double)(nseg - 1);
if (smooth_option == SMOOTH_ACSPLINES)
sc_z = cp_approx_spline(curve->points, curve->p_count, 0, 2, 3);
else
sc_z = cp_tridiag(curve->points, curve->p_count, 0, 2);
for (i = 0, l = 0; i < nseg; i++) {
double temp;
t = old_points[0].x + i * tstep;
/* Move forward to the spline interval this point is in */
while ((t >= old_points[l + 1].x) && (l < curve->p_count- 2))
l++;
temp = t - old_points[l].x;
new_points[i].x = t;
new_points[i].y = old_points[l].y; /* All the same */
new_points[i].z = ((sc_z[l][3] * temp + sc_z[l][2]) * temp + sc_z[l][1])
* temp + sc_z[l][0];
}
}
/*
* Check for a curve lying in the xy plane (z = constant).
*/
else if (maxdz < FLT_EPSILON) {
tstep = (old_points[curve->p_count-1].x - old_points[0].x) / (double)(nseg - 1);
if (smooth_option == SMOOTH_ACSPLINES)
sc_y = cp_approx_spline(curve->points, curve->p_count, 0, 1, 3);
else
sc_y = cp_tridiag(curve->points, curve->p_count, 0, 1);
for (i = 0, l = 0; i < nseg; i++) {
double temp;
t = old_points[0].x + i * tstep;
/* Move forward to the spline interval this point is in */
while ((t >= old_points[l + 1].x) && (l < curve->p_count- 2))
l++;
temp = t - old_points[l].x;
new_points[i].x = t;
new_points[i].y = ((sc_y[l][3] * temp + sc_y[l][2]) * temp + sc_y[l][1])
* temp + sc_y[l][0];
new_points[i].z = old_points[l].z; /* All the same */
}
}
/*
* This is the general case.
* Calculate spline coefficients for each dimension x, y, z
*/
else {
if (smooth_option == SMOOTH_ACSPLINES) {
sc_x = cp_approx_spline(curve->points, curve->p_count, PATHCOORD, 0, 3);
sc_y = cp_approx_spline(curve->points, curve->p_count, PATHCOORD, 1, 3);
sc_z = cp_approx_spline(curve->points, curve->p_count, PATHCOORD, 2, 3);
} else {
sc_x = cp_tridiag( curve->points, curve->p_count, PATHCOORD, 0);
sc_y = cp_tridiag( curve->points, curve->p_count, PATHCOORD, 1);
sc_z = cp_tridiag( curve->points, curve->p_count, PATHCOORD, 2);
}
for (i = 0, l=0; i < nseg; i++) {
double temp;
t = i * tstep;
/* Move forward to the spline interval this point is in */
while ((t >= old_points[l + 1].CRD_PATH) && (l < curve->p_count- 2))
l++;
temp = t - old_points[l].CRD_PATH;
new_points[i].x = ((sc_x[l][3] * temp + sc_x[l][2]) * temp + sc_x[l][1])
* temp + sc_x[l][0];
new_points[i].y = ((sc_y[l][3] * temp + sc_y[l][2]) * temp + sc_y[l][1])
* temp + sc_y[l][0];
new_points[i].z = ((sc_z[l][3] * temp + sc_z[l][2]) * temp + sc_z[l][1])
* temp + sc_z[l][0];
}
}
/* We're done with the spline coefficients */
free(sc_x);
free(sc_y);
free(sc_z);
/* Replace original data with spline approximation */
free(curve->points);
curve->points = new_points;
curve->p_count = nseg;
curve->p_max = nseg+1; /* not sure why we asked for 1 extra */
}
/*
* Generate 2D splines along a path for each set of points in the plot,
* smoothing option SMOOTH_PATH.
* TODO:
* - number of spline samples should be controlled by something other
* than "set samples"
* - spline weights from an additional column
*/
void
gen_2d_path_splines( struct curve_points *plot )
{
int i;
int ic, first_point; /* indexes for original data */
int is = 0; /* index for new (splined) data */
struct coordinate *old_points = NULL;
struct coordinate *splined_points;
spline_coeff *sc_x = NULL;
spline_coeff *sc_y = NULL;
double xrange = fabs(axis_array[plot->x_axis].max - axis_array[plot->x_axis].min);
double yrange = fabs(axis_array[plot->y_axis].max - axis_array[plot->y_axis].min);
int curves = num_curves(plot);
/* Allocate space to hold the interpolated points */
splined_points = gp_alloc( (plot->p_count + samples_1 * curves) * sizeof(struct coordinate), NULL );
memset( splined_points, 0, (plot->p_count + samples_1 * curves) * sizeof(struct coordinate));
first_point = 0;
for (ic = 0; ic < curves; ic++) {
double t, tstep, tsum;
double dx, dy;
int l;
int nold;
int num_points = next_curve(plot, &first_point);
TBOOLEAN closed = FALSE;
/* Make a copy of the original points so that we don't corrupt the
* list by adding up to three new ones.
*/
old_points = gp_realloc( old_points, (num_points + 3) * sizeof(struct coordinate),
"spline points");
memcpy( &old_points[1], &plot->points[first_point], num_points * sizeof(struct coordinate));
/* Remove any unusable points (NaN, missing, duplicates) before fitting a spline.
* If that leaves fewer than 3 points, skip it.
*/
nold = do_curve_cleanup(&old_points[1], num_points);
if (nold < 3) {
first_point += num_points;
continue;
}
/* We expect one of two cases. Either this really is a closed
* curve (end point matches start point) or it is an open-ended
* path that may not be monotonic on x.
* For plot style "with filledcurves closed" we add an extra
* point at the end if it is not already there.
*/
if (old_points[1].x == old_points[nold].x
&& old_points[1].y == old_points[nold].y)
closed = TRUE;
if ((plot->plot_style == FILLEDCURVES) && !closed) {
old_points[++nold] = old_points[1];
closed = TRUE;
}
if (closed) {
/* Wrap around to one point before and one point after the path closure */
nold += 2;
old_points[0] = old_points[nold-3];
old_points[nold-1] = old_points[2];
} else {
/* Dummy up an extension at either end */
nold += 2;
old_points[0].x = old_points[1].x + old_points[1].x - old_points[2].x;
old_points[nold-1].x = old_points[nold-2].x + old_points[nold-2].x - old_points[nold-3].x;
old_points[0].y = old_points[1].y + old_points[1].y - old_points[2].y;
old_points[nold-1].y = old_points[nold-2].y + old_points[nold-2].y - old_points[nold-3].y;
}
/* Construct path-length vector; store it in an unused slot of old_points */
t = tsum = 0.0;
old_points[0].CRD_PATH = 0;
if (xrange == 0)
xrange = 1.;
if (yrange == 0)
yrange = 1.;
for (i = 1; i < nold; i++) {
dx = (old_points[i].x - old_points[i-1].x) / xrange;
dy = (old_points[i].y - old_points[i-1].y) / yrange;
tsum += sqrt( dx*dx + dy*dy );
old_points[i].CRD_PATH = tsum;
}
/* Normalize so that the path fraction always runs from 0 to 1 */
for (i = 1; i < nold; i++)
old_points[i].CRD_PATH /= tsum;
tstep = 1.0 / (double)(samples_1 - 1);
/* Calculate spline coefficients for x and for y as a function of path */
sc_x = cp_tridiag( old_points, nold, PATHCOORD, 0);
sc_y = cp_tridiag( old_points, nold, PATHCOORD, 1);
/* First output point is the same as the original first point */
splined_points[is++] = old_points[1];
/* Skip the points in the overlap region */
for (i = 0; i * tstep < old_points[1].CRD_PATH; i++)
;
/* Use spline coefficients to generate a new point at each sample interval. */
for (l=0; i < samples_1; i++) {
double temp;
t = i * tstep;
/* Stop before wrapping around. Copy the original end point. */
if (t > old_points[nold-2].CRD_PATH) {
splined_points[is++] = old_points[nold-2];
break;
}
/* Move forward to the spline interval this point is in */
while ((t >= old_points[l + 1].CRD_PATH) && (l < nold- 2))
l++;
temp = t - old_points[l].CRD_PATH;
splined_points[is].x = ((sc_x[l][3] * temp + sc_x[l][2]) * temp + sc_x[l][1])
* temp + sc_x[l][0];
splined_points[is].y = ((sc_y[l][3] * temp + sc_y[l][2]) * temp + sc_y[l][1])
* temp + sc_y[l][0];
is++;
}
/* Done with spline coefficients */
free(sc_x);
free(sc_y);
/* Add a seperator point after this set of splined points */
splined_points[is++].type = UNDEFINED;
first_point += num_points;
}
/* Replace original data with splined approximation */
free(old_points);
free(plot->points);
plot->points = splined_points;
plot->p_max = curves * samples_1;
plot->p_count = is;
return;
}
/*
* Externally callable interface to 3D spline routines
*/
void
gen_3d_splines( struct surface_points *plot )
{
struct iso_curve *curve = plot->iso_crvs;
while (curve) {
/* Remove any unusable points before fitting a spline */
curve->p_count = do_curve_cleanup(curve->points, curve->p_count);
if (curve->p_count > 3)
do_3d_cubic(curve, plot->plot_smooth);
curve = curve->next;
}
}
static int
do_curve_cleanup( struct coordinate *point, int npoints )
{
int i, keep;
/* Step through points in curve keeping only the usable ones.
* Discard duplicates
*/
keep = 0;
for (i = 0; i < npoints; i++) {
if (point[i].type == UNDEFINED)
continue;
if (isnan(point[i].x) || isnan(point[i].y) || isnan(point[i].z))
continue;
if (i != keep)
point[keep] = point[i];
/* FIXME: should probably check fabs(this-prev) < EPS */
if ((keep > 0) && (point[keep].x == point[keep-1].x)
&& (point[keep].y == point[keep-1].y)
&& (point[keep].z == point[keep-1].z))
continue;
keep++;
}
return keep;
}
/*
* convex_hull() replaces the original set of points with a subset that
* delimits the convex hull of the original points.
* The hull is found using Graham's algorithm
* RL Graham (1972), Information Processing Letters 1: 132133.
* winnow_interior_points() is a helper routine that can greatly reduce
* processing time for large data sets but is otherwise not necessary.
* expand_hull() pushes the hull segments away from the interior to
* produce a bounding curve exterior to all points.
* - Ethan A Merritt 2021
*/
#define CROSS(p1,p2,p3) \
( ((p2)->x - (p1)->x) * ((p3)->y - (p2)->y) \
- ((p2)->y - (p1)->y) * ((p3)->x - (p2)->x) )
void
convex_hull(struct curve_points *plot)
{
int i;
struct coordinate *points = plot->points;
struct coordinate **stack = NULL;
int np = plot->p_count;
int ntop;
/* cluster is used to hold x/y ranges and centroid */
t_cluster cluster;
cluster.npoints = plot->p_count;
/* Special cases */
if (np < 3)
return;
if (np == 3) {
cp_extend(plot, 4);
plot->points[3] = plot->points[0];
plot->p_count = 4;
return;
}
/* Find x and y limits of points in the cluster. */
cluster_stats(plot->points, &cluster);
/* This is not strictly necessary, but greatly reduces the number
* of points to be sorted and tested for the hull boundary.
*/
winnow_interior_points(plot, &cluster);
/* Sort the remaining points (probably only need to sort on x?) */
gp_qsort(plot->points, plot->p_count,
sizeof(struct coordinate), compare_xyz);
/* Find hull points using a variant of Graham's algorithm.
* The path through the points is accumulated on a stack.
*/
stack = gp_alloc( (np+1)*sizeof(void *), "Hull" );
/* Initialize stack with known start of top arc and first candidate point */
stack[0] = &points[0];
stack[1] = &points[1];
np = 2;
for (i=2; ip_count; i++) {
while ((np >= 2) && CROSS( stack[np-2], stack[np-1], &points[i] ) >= 0)
np--;
stack[np++] = &points[i];
}
/* push onto stack the first candidate point for lower arc */
i -= 2;
stack[np++] = &points[i];
ntop = np;
for (i--; i>=0; i--) {
while ( (np >= ntop) && CROSS( stack[np-2], stack[np-1], &points[i] ) >= 0)
np--;
stack[np++] = &points[i];
}
/* Replace the original list of points with the ordered path */
points = gp_alloc( np * sizeof(struct coordinate), "Hull" );
for (i=0; ipoints = points;
plot->p_count = np;
plot->p_max = np;
}
/*
* winnow_interior_points() is an optional helper routine for convex_hull.
* It reduces the number of points to be sorted and processed by removing
* points in a quadrilateral bounded by the four points with max/min x/y.
*/
static void
winnow_interior_points (struct curve_points *plot, t_cluster *cluster)
{
#define TOLERANCE -1.e-10
struct coordinate *p, *pp1, *pp2, *pp3, *pp4;
struct coordinate *points = plot->points;
double area;
int i, np;
/* Find points defining maximal extent on x and y */
pp1 = pp2 = pp3 = pp4 = plot->points;
for (p = plot->points; p < &(plot->points[plot->p_count]); p++) {
if (p->x == cluster->xmin) pp1 = p;
if (p->x == cluster->xmax) pp2 = p;
if (p->y == cluster->ymin) pp3 = p;
if (p->y == cluster->ymax) pp4 = p;
}
/* Ignore any points that lie inside the clockwise triangle bounded by pp1 pp2 pp3 */
area = fabs(-pp2->y*pp3->x + pp1->y*(-pp2->x + pp3->x)
+ pp1->x*(pp2->y - pp3->y) + pp2->x*pp3->y);
area += TOLERANCE;
for (i=0; ip_count; i++) {
double px = points[i].x;
double py = points[i].y;
double s = (pp1->y*pp3->x - pp1->x*pp3->y + (pp3->y - pp1->y)*px + (pp1->x - pp3->x)*py);
double t = (pp1->x*pp2->y - pp1->y*pp2->x + (pp1->y - pp2->y)*px + (pp2->x - pp1->x)*py);
if ( (s < TOLERANCE) && (t < TOLERANCE) && (fabs(s+t) < area) )
points[i].type = EXCLUDEDRANGE;
}
/* Also ignore points in the clockwise triangle bounded by pp3 pp4 pp1 */
area = fabs(-pp4->y*pp1->x + pp3->y*(-pp4->x + pp1->x)
+ pp3->x*(pp4->y - pp1->y) + pp4->x*pp1->y);
area += TOLERANCE;
for (i=0; ip_count; i++) {
double px = points[i].x;
double py = points[i].y;
double s = (pp3->y*pp1->x - pp3->x*pp1->y + (pp1->y - pp3->y)*px + (pp3->x - pp1->x)*py);
double t = (pp3->x*pp4->y - pp3->y*pp4->x + (pp3->y - pp4->y)*px + (pp4->x - pp3->x)*py);
if ( (s < TOLERANCE) && (t < TOLERANCE) && (fabs(s+t) < area) )
points[i].type = EXCLUDEDRANGE;
}
/* Discard the interior points and outliers */
np = 0;
p = points;
for (i=0; ip_count; i++) {
if (points[i].type == UNDEFINED)
continue;
if (points[i].type != EXCLUDEDRANGE)
p[np++] = points[i];
}
plot->p_count = np;
#undef TOLERANCE
}
/*
* expand_hull() "inflates" a convex or concave hull by displacing
* each edge away from the interior along its normal vector
* by requested distance d.
*/
void
expand_hull(struct curve_points *plot)
{
struct coordinate *newpoints = NULL;
struct coordinate *points = plot->points;
double scale = plot->smooth_parameter;
double d = fabs(scale);
int N = plot->p_count;
struct coordinate *v1, *v2, *v3;
t_cluster cluster;
double winding; /* +1 for clockwise -1 for anticlockwise */
int ni;
/* Determine whether the hull points are ordered clockwise or
* anticlockwise. This allows easy determination of the interior
* side of each edge and the concave/convex status of each vertex.
*/
cluster.npoints = N;
cluster_stats(points, &cluster);
v1 = &points[ (cluster.pin == 0) ? N-2 : cluster.pin-1 ];
v2 = &points[cluster.pin];
v3 = &points[cluster.pin + 1];
if (CROSS(v2,v1,v3) > 0)
winding = 1.0;
else
winding = -1.0;
/* Each edge of the hull is displaced outward by a constant amount d.
* At each convex vertex this replaces the original point with two
* points that are the vertices of a beveled join.
* At each concave vertex this displaces the original point toward
* the mouth of the concavity.
*/
newpoints = gp_alloc(2*N * sizeof(struct coordinate), "expand hull");
ni = 0;
for (int i = 0; i < N; i++) {
double m1, m2, d2norm;
double dx1, dy1, dx2, dy2;
v1 = &points[ (i == 0) ? N-2 : i-1 ];
v2 = &points[i];
v3 = &points[ (i == N-1) ? 1 : i+1 ];
m1 = (v2->y - v1->y) / (v2->x - v1->x); /* slope of v1v2 */
d2norm = d*d / (1/(m1*m1) + 1);
dx1 = copysign( sqrt(d2norm), -winding * (v2->y - v1->y) );
dy1 = copysign( sqrt(d*d - d2norm), winding * (v2->x - v1->x) );
m2 = (v3->y - v2->y) / (v3->x - v2->x); /* slope of v2v3 */
d2norm = d*d / (1/(m2*m2) + 1);
dx2 = copysign( sqrt(d2norm), -winding * (v3->y - v2->y) );
dy2 = copysign( sqrt(d*d - d2norm), winding * (v3->x - v2->x) );
/* convex vertex */
if (winding * CROSS(v1,v2,v3) < 0) {
newpoints[ni] = points[i];
newpoints[ni].x = v2->x + dx1;
newpoints[ni].y = v2->y + dy1;
ni++;
newpoints[ni] = points[i];
newpoints[ni].x = v2->x + dx2;
newpoints[ni].y = v2->y + dy2;
ni++;
/* concave vertex (over-emphasizes steep holes) */
} else {
double dnorm = d / sqrt((dx1+dx2)*(dx1+dx2) + (dy1+dy2)*(dy1+dy2));
newpoints[ni] = points[i];
newpoints[ni].x = v2->x + dnorm * (dx1 + dx2);
newpoints[ni].y = v2->y + dnorm * (dy1 + dy2);
ni++;
}
}
/* Replace original point list with the new one */
cp_extend(plot, 0);
plot->points = newpoints;
plot->p_count = ni;
plot->p_max = 2*N;
}
/*
* Find min, max, center of mass points in cluster held in plot->points.
* All points not marked EXCLUDEDRANGE are assumed to be in a single cluster.
*/
static void
cluster_stats( struct coordinate *points, t_cluster *cluster )
{
double xsum = 0;
double ysum = 0;
int excluded = 0;
cluster->xmin = cluster->ymin = VERYLARGE;
cluster->xmax = cluster->ymax = -VERYLARGE;
cluster->pin = -1;
for (int i = 0; i < cluster->npoints; i++) {
if (points[i].type == EXCLUDEDRANGE || points[i].type == UNDEFINED) {
excluded++;
continue;
}
xsum += points[i].x;
ysum += points[i].y;
if ((cluster->xmin == points[i].x)
&& (points[i].y > points[cluster->pin].y))
cluster->pin = i;
if (cluster->xmin > points[i].x) {
cluster->xmin = points[i].x;
cluster->pin = i;
}
if (cluster->ymin > points[i].y)
cluster->ymin = points[i].y;
if (cluster->xmax < points[i].x)
cluster->xmax = points[i].x;
if (cluster->ymax < points[i].y)
cluster->ymax = points[i].y;
}
cluster->cx = xsum / (cluster->npoints - excluded);
cluster->cy = ysum / (cluster->npoints - excluded);
}
/*
* The "sharpen" filter looks for truncated extrema in the function being plotted.
* The true local extremum is found by bisection and added to the set of
* points being plotted.
*
* Motivation:
* If the function being plotted has a sharp extremum that lies between
* two of the sampled points, the resulting plot truncates the extremum.
* Increasing the number of samples would help, of course, but they are
* only really needed close to the extremum. And no matter how many
* samples there are, the true peak may lie between two of them.
* Example: plot lgamma(x)
*
* Method:
* - Approximate the second derivative of the function at each point using
* a 5-term Savitsky-Golay filter.
* - Inspect a moving window centered on point i
* Truncation is suspected if either
* f(i) is a minimum and f''(i-2) > 0 and/or f''(i+2) > 0
* or
* f(i) is a maximum and f''(i-2) < 0 and/or f''(i+2) < 0
* - Use bisection to find the true local extremum
* - Add the new point next to the original point i
*
* Ethan A Merritt Dec 2022
*
* TODO:
* - Diagnose and warn that finer sampling might help?
*/
void
sharpen(struct curve_points *plot)
{
struct axis *y_axis = &axis_array[plot->y_axis];
int newcount = plot->p_count;
struct coordinate *p;
/* Restrictions */
if (parametric || polar)
return;
/* Make more than enough room for new points */
cp_extend(plot, 1.5 * plot->p_count);
p = plot->points;
/* We expect that sharpened peaks may go to +/- Infinity
* and in fact the plot may already contain such points.
* Because these were tagged UNDEFINED in store_and_update_range()
* they will be lost when the points are sorted.
* As a work-around we look for these and flag them OUTRANGE instead.
* If there are other UNDEFINED points we set y to 0 so that they
* do not cause the f'' approximation to blow up.
*/
for (int i = 0; i < plot->p_count; i++) {
if (p[i].type == UNDEFINED) {
if (p[i].y >= VERYLARGE) {
p[i].type = OUTRANGE;
p[i].y = VERYLARGE;
} else if (p[i].y p_count-4; i++) {
TBOOLEAN criterion = FALSE;
criterion = (fpp_SG5(&p[i-2]) < 0 || fpp_SG5(&p[i+2]) < 0);
if (p[i].y varcolor)
for (i = 0; i < plot->p_count; i++)
plot->points[i].CRD_COLOR = plot->varcolor[i];
/* Sort the points on x.
* This allows us to reduce the time required from O(N^2) to
* approximately O(NlogN).
*/
gp_qsort(plot->points, N, sizeof(struct coordinate), compare_xyz);
/* Construct a triangle "sufficiently big" to enclose the set of points.
* That means each bounding vertex must be far enough away from the data
* points that a circumcircle does not catch extra data points.
* It turns out this is easier to guarantee if we start with two triangles.
*/
{
double xmin = VERYLARGE, xmax = -VERYLARGE;
double ymin = VERYLARGE, ymax = -VERYLARGE;
double xdelta, ydelta, gap;
/* Allocate 4 additional points for the two triangles */
cp_extend(plot, N+4);
points = plot->points;
for (p = points; p < &points[N]; p++) {
if (p->type == UNDEFINED)
continue;
if (p->x < xmin) xmin = p->x;
if (p->x > xmax) xmax = p->x;
if (p->y < ymin) ymin = p->y;
if (p->y > ymax) ymax = p->y;
}
xdelta = xmax - xmin;
ydelta = ymax - ymin;
/* Is 100 "sufficiently big"??? */
gap = 100. * GPMAX(xdelta,ydelta);
points[N].x = xmin - gap;
points[N+1].x = xmax + gap;
points[N+2].x = xmin - gap;
points[N+3].x = xmax + gap;
points[N].y = ymin - gap;
points[N+1].y = ymin - gap;
points[N+2].y = ymax + gap;
points[N+3].y = ymax + gap;
t = insert_new_triangle( points, N, N+1, N+2 );
find_circumcircle(points, t);
t = insert_new_triangle( points, N+1, N+2, N+3 );
find_circumcircle(points, t);
}
/* Add points one by one */
for (p = points; p < &points[N]; p++) {
/* Ignore undefined points. */
if (p->type == UNDEFINED)
continue;
/* Also ignore duplicate points */
if (p->x == (p+1)->x && p->y == (p+1)->y)
continue;
/* First step is to move all triangles for which the new point
* violates the criterion "no other points in bounding circle"
* to a separate list bad_triangles.
*/
delete_triangles(&bad_triangles);
prev = &good_triangles;
for (t = prev->next; t; t = prev->next) {
if (in_circumcircle(t, p->x, p->y))
invalidate_triangle(prev, t);
/* Points are sorted on x, so if the circumcircle of this
* triangle is entirely to the left of the current point
* we needn't check it for subsequent points either.
*/
else if ((t->cx + t->r) < p->x) {
freeze_triangle(prev, t);
if (frozen_tail == NULL)
frozen_tail = t;
} else
prev = t;
}
/* I think this only happens if there is a duplicate point */
if (bad_triangles.next == NULL)
continue;
/* Second step is to find the edges forming the perimeter
* of the bad triangles. An edge cannot appear twice.
* If it does, flag both copies as -1.
*/
for (i = 0, t = bad_triangles.next; t; t = t->next)
i++;
edges = gp_realloc(edges, 3*i * sizeof(t_edge), "delaunay edges");
for (i = 0, t = bad_triangles.next; t; t = t->next) {
edges[i].v1 = (t->v1 > t->v2) ? t->v2 : t->v1;
edges[i].v2 = (t->v1 > t->v2) ? t->v1 : t->v2;
i++;
edges[i].v1 = (t->v2 > t->v3) ? t->v3 : t->v2;
edges[i].v2 = (t->v2 > t->v3) ? t->v2 : t->v3;
i++;
edges[i].v1 = (t->v3 > t->v1) ? t->v1 : t->v3;
edges[i].v2 = (t->v3 > t->v1) ? t->v3 : t->v1;
i++;
}
nedges = i;
gp_qsort(edges, nedges, sizeof(t_edge), compare_edges);
for (i = 0; i < nedges-1; i++) {
if (edges[i].v1 == edges[i+1].v1 && edges[i].v2 == edges[i+1].v2)
edges[i].v1 = edges[i+1].v1 = -1;
}
/* For each edge on the perimeter construct a new triangle
* containing the new point as its third vertex.
*/
for (i = 0; i < nedges; i++) {
if (edges[i].v1 >= 0) {
t = insert_new_triangle( points, edges[i].v1, edges[i].v2, p - points);
find_circumcircle(points, t);
}
}
}
/* Merge frozen triangles back into the good triangle list */
if (frozen_tail) {
frozen_tail->next = good_triangles.next;
good_triangles.next = frozen_triangles.next;
frozen_triangles.next = NULL;
frozen_tail = NULL;
}
/* Remove any triangles that contain a vertex of the original "big triangle".
* The perimeter of the resulting tesselation can be found later by
* collecting one edge from each of the triangles left in the bad_triangles list.
*/
delete_triangles(&bad_triangles);
prev = &good_triangles;
for (t = prev->next; t; t = prev->next) {
if ( (t->v1 >= N) || (t->v2 >= N) || (t->v3 >= N) )
invalidate_triangle(prev, t);
else
prev = t;
}
/* Flag the original points that lie on the bounding hull.
* The edges making up the hull remain in bounding_edges[]
* for use by a subsequent call to concave_hull().
*/
for (i = 0; i < N; i++)
points[i].extra = 0;
for (i = 0, t = bad_triangles.next; t; t = t->next)
i++;
edges = gp_realloc(edges, 3*i * sizeof(t_edge), "delaunay edges");
max_bounding_edges = 3*i;
nedges = 0;
for (t = bad_triangles.next; t; t = t->next) {
int v1, v2;
v1 = (t->v1 >= N) ? t->v3 : t->v1;
v2 = (t->v2 >= N) ? t->v3 : t->v2;
if ((v1 == v2) || (v1 >= N) || (v2 >= N))
continue;
edges[nedges].v1 = v1;
edges[nedges].v2 = v2;
points[v1].extra = HULL_POINT;
points[v2].extra = HULL_POINT;
edges[nedges].length = edge_length( &edges[nedges], plot->points );
nedges++;
}
bounding_edges = edges;
n_bounding_edges = nedges;
}
static void
invalidate_triangle( triangle *prev, triangle *this )
{
prev->next = this->next;
this->next = bad_triangles.next;
bad_triangles.next = this;
}
static void
freeze_triangle( triangle *prev, triangle *this )
{
prev->next = this->next;
this->next = frozen_triangles.next;
frozen_triangles.next = this;
}
static int
compare_edges(SORTFUNC_ARGS e1, SORTFUNC_ARGS e2)
{
t_edge *t1 = (t_edge *)e1;
t_edge *t2 = (t_edge *)e2;
if (t1->v1 > t2->v1)
return 1;
if (t1->v1 < t2->v1)
return -1;
if (t1->v2 > t2->v2)
return 1;
if (t1->v2 < t2->v2)
return -1;
return 0;
}
static triangle *
insert_new_triangle( struct coordinate *points, int v1, int v2, int v3 )
{
triangle *t = gp_alloc( sizeof(triangle), "triangle" );
if (v1 < v2 && v1 < v3) {
t->v1 = v1;
t->v2 = (v2 < v3) ? v2 : v3;
t->v3 = (v2 < v3) ? v3 : v2;
} else if (v2 < v1 && v2 < v3) {
t->v1 = v2;
t->v2 = (v1 < v3) ? v1 : v3;
t->v3 = (v1 < v3) ? v3 : v1;
} else {
t->v1 = v3;
t->v2 = (v1 < v2) ? v1 : v2;
t->v3 = (v1 < v2) ? v2 : v1;
}
t->next = good_triangles.next;
good_triangles.next = t;
return t;
}
static void
find_circumcircle( struct coordinate *points, triangle *t)
{
double BB_x, BB_y, CC_x, CC_y, DD;
double cent_x, cent_y;
BB_x = points[t->v2].x - points[t->v1].x;
BB_y = points[t->v2].y - points[t->v1].y;
CC_x = points[t->v3].x - points[t->v1].x;
CC_y = points[t->v3].y - points[t->v1].y;
DD = 2. * (BB_x * CC_y - BB_y * CC_x);
cent_x = CC_y * (BB_x*BB_x + BB_y*BB_y) - BB_y * (CC_x*CC_x + CC_y*CC_y);
cent_x /= DD;
cent_y = BB_x * (CC_x*CC_x + CC_y*CC_y) - CC_x * (BB_x*BB_x + BB_y*BB_y);
cent_y /= DD;
t->r = sqrt( cent_x * cent_x + cent_y * cent_y );
t->cx = cent_x + points[t->v1].x;
t->cy = cent_y + points[t->v1].y;
}
static TBOOLEAN
in_circumcircle( triangle *t, double x, double y )
{
if ((t->cx - x)*(t->cx - x) + (t->cy - y)*(t->cy - y) < (t->r * t->r))
return TRUE;
else
return FALSE;
}
static void
delete_triangles( triangle *list_head )
{
triangle *next;
triangle *t = list_head->next;
while (t) {
next = t->next;
free(t);
t = next;
}
list_head->next = NULL;
}
static int
compare_edgelength(SORTFUNC_ARGS e1, SORTFUNC_ARGS e2)
{
if (((t_edge *)e1)->length > ((t_edge *)e2)->length)
return -1;
if (((t_edge *)e1)->length < ((t_edge *)e2)->length)
return 1;
return 0;
}
/* Adjust current perimeter of triangulated points by incrementally
* removing single edges to reveal a portion of the interior.
* Assume
* plot->points contains the data points
* bounding_edges[] contains the convex hull
* l_chi is the characteristic length determining how far we trim
* Iterative algorithm
* Sort bounding edges by length, i.e. bounding_edges[0] is the longest
* If there is no removable edge with length > l_chi, stop.
* Check whether longest edge belongs to a triangle with 3 exterior points
* If so
* mark edge as non-removable and continue
* If not
* Remove it to expose two new edges
* Add the two newly exposed edges to the perimeter.
* Remove triangle from good_triangles list.
* Note
* On exit some lengths in bounding_edges have been overwritten with -1
*/
static void
chi_reduction( struct curve_points *plot, double l_chi )
{
/* Each edge removed exposes two others. We will increase max as necessary. */
bounding_edges = gp_realloc( bounding_edges, 2*n_bounding_edges*sizeof(t_edge),
"bounding_edges");
max_bounding_edges = 2*n_bounding_edges;
gp_qsort(bounding_edges, n_bounding_edges, sizeof(t_edge), compare_edgelength);
while (dig_one_hole( plot, l_chi )) {
gp_qsort(bounding_edges, n_bounding_edges, sizeof(t_edge), compare_edgelength);
}
}
static double
edge_length( t_edge *edge, struct coordinate *p )
{
double dx = p[edge->v1].x - p[edge->v2].x;
double dy = p[edge->v1].y - p[edge->v2].y;
return sqrt(dx*dx + dy*dy);
}
static void
keep_edge( int i, int v1, int v2, struct coordinate *points)
{
bounding_edges[i].v1 = v1;
bounding_edges[i].v2 = v2;
bounding_edges[i].length = edge_length(&bounding_edges[i], points);
}
static TBOOLEAN
dig_one_hole( struct curve_points *plot, double l_chi )
{
struct coordinate *p = plot->points;
triangle *t, *prev;
if (bounding_edges[0].length next; t; t = t->next) {
TBOOLEAN hit = FALSE;
if (bounding_edges[0].v1 == t->v1) {
if (bounding_edges[0].v2 == t->v2) {
hit = TRUE;
if (p[t->v3].extra == HULL_POINT) {
bounding_edges[0].length = -1; /* Flag edge as not removeable */
} else {
p[t->v3].extra = HULL_POINT; /* Mark uncovered point exterior */
keep_edge(0, t->v1, t->v3, p);
keep_edge(n_bounding_edges, t->v2, t->v3, p);
n_bounding_edges++;
invalidate_triangle(prev, t);
}
} else if (bounding_edges[0].v2 == t->v3) {
hit = TRUE;
if (p[t->v2].extra == HULL_POINT) {
bounding_edges[0].length = -1; /* Flag edge as not removeable */
} else {
p[t->v2].extra = HULL_POINT; /* Mark uncovered point exterior */
keep_edge(0, bounding_edges[0].v1, t->v2, p);
keep_edge( n_bounding_edges, t->v2, t->v3, p);
n_bounding_edges++;
invalidate_triangle(prev, t);
}
}
} else if (bounding_edges[0].v1 == t->v2) {
if (bounding_edges[0].v2 == t->v3) {
hit = TRUE;
if (p[t->v1].extra == HULL_POINT) {
bounding_edges[0].length = -1; /* Flag edge as not removeable */
} else {
p[t->v1].extra = HULL_POINT; /* Mark uncovered point exterior */
keep_edge(0, t->v1, t->v2, p);
keep_edge(n_bounding_edges, t->v1, t->v3, p);
n_bounding_edges++;
invalidate_triangle(prev, t);
}
}
}
/* Replaced old perimeter segment with two new ones */
if (hit) {
if (n_bounding_edges >= max_bounding_edges) {
max_bounding_edges *= 2;
bounding_edges = gp_realloc(bounding_edges, max_bounding_edges * sizeof(t_edge),
"bounding_edges");
}
return(TRUE);
}
/* Otherwise continue to search for triangle containing this edge */
prev = t;
}
return(FALSE);
}
void
concave_hull( struct curve_points *plot )
{
udvt_entry *udv;
int prev;
struct coordinate *newpoints;
double chi_length = 0.0;
/* Construct -shape */
udv = get_udv_by_name("chi_length");
if (udv && udv->udv_value.type == CMPLX)
chi_length = real(&(udv->udv_value));
/* No chi_length given, make a guesstimate */
if (chi_length points[ bounding_edges[0].v1 ];
newpoints[1] = plot->points[ bounding_edges[0].v2 ];
prev = bounding_edges[0].v2;
bounding_edges[0].v1 = bounding_edges[0].v2 = -1;
for (int n = 2; n < n_bounding_edges; n++) {
for (int i = 1; i < n_bounding_edges; i++) {
if (prev == bounding_edges[i].v1)
prev = bounding_edges[i].v2;
else if (prev == bounding_edges[i].v2)
prev = bounding_edges[i].v1;
else /* Not a match; keep looking */
continue;
newpoints[n] = plot->points[prev];
bounding_edges[i].v1 = bounding_edges[i].v2 = -1;
break;
}
}
newpoints[n_bounding_edges] = newpoints[0];
cp_extend(plot, 0);
plot->p_max = n_bounding_edges+1;
plot->p_count = n_bounding_edges+1;
plot->points = newpoints;
}
/*
* filter option "delaunay"
* plot POINTS using 1:2 delaunay with polygons fs empty
* Replace original list of points with a polygon representation
* of the individual triangles.
*
* This routine was written to help debug implementation of -shapes
* to generate a concave hull. See "concave_hull.dem".
* It is currently undocumented as a separate filter option.
*/
void
save_delaunay_triangles( struct curve_points *plot )
{
struct coordinate *point = plot->points;
struct coordinate *newpoints = NULL;
double *newcolor = NULL;
triangle *t;
int outp;
/* Reserve space to store each triangle as three points
* plus one "empty" spacer.
*/
for (outp = 0, t = good_triangles.next; t; t = t->next)
outp++;
newpoints = gp_alloc( 5 * outp * sizeof(struct coordinate), "delaunay filter" );
/* Copy full original point so that it retains any extra fields */
for (outp = 0, t = good_triangles.next; t; t = t->next) {
newpoints[outp++] = point[t->v1];
newpoints[outp++] = point[t->v2];
newpoints[outp++] = point[t->v3];
newpoints[outp++] = point[t->v1];
newpoints[outp++] = blank_data_line;
}
/* Make a new separate list of colors */
if (plot->varcolor) {
newcolor = gp_alloc( 5 * outp * sizeof(double), "delaunay colors" );
for (int i = 0; i < outp; i++)
newcolor[i] = newpoints[i].CRD_COLOR;
}
/* Replace original point list with the new one */
cp_extend(plot, 0);
plot->points = newpoints;
plot->varcolor = newcolor;
plot->p_count = outp;
plot->p_max = outp;
}
#endif /* WITH_CHI_SHAPES */