1/*
2 * "$Id: curve.c,v 1.55 2010/08/04 00:33:56 rlk Exp $"
3 *
4 *   Print plug-in driver utility functions for the GIMP.
5 *
6 *   Copyright 1997-2000 Michael Sweet (mike@easysw.com) and
7 *	Robert Krawitz (rlk@alum.mit.edu)
8 *
9 *   This program is free software; you can redistribute it and/or modify it
10 *   under the terms of the GNU General Public License as published by the Free
11 *   Software Foundation; either version 2 of the License, or (at your option)
12 *   any later version.
13 *
14 *   This program is distributed in the hope that it will be useful, but
15 *   WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 *   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
17 *   for more details.
18 *
19 *   You should have received a copy of the GNU General Public License
20 *   along with this program; if not, write to the Free Software
21 *   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22 */
23
24#ifdef HAVE_CONFIG_H
25#include <config.h>
26#endif
27#include <gutenprint/gutenprint.h>
28#include "gutenprint-internal.h"
29#include <gutenprint/gutenprint-intl-internal.h>
30#include <math.h>
31#ifdef sun
32#include <ieeefp.h>
33#endif
34#include <string.h>
35#include <stdlib.h>
36#include <limits.h>
37#include <unistd.h>
38#include <sys/types.h>
39#include <sys/stat.h>
40
41#ifdef __GNUC__
42#define inline __inline__
43#endif
44
45#undef inline
46#define inline
47
48/*
49 * We could do more sanity checks here if we want.
50 */
51#define CHECK_CURVE(curve)			\
52do						\
53  {						\
54    STPI_ASSERT((curve) != NULL, NULL);		\
55    STPI_ASSERT((curve)->seq != NULL, NULL);	\
56  } while (0)
57
58static const int curve_point_limit = 1048576;
59
60struct stp_curve
61{
62  stp_curve_type_t curve_type;
63  stp_curve_wrap_mode_t wrap_mode;
64  int piecewise;
65  int recompute_interval;	/* Do we need to recompute the deltas? */
66  double gamma;			/* 0.0 means that the curve is not a gamma */
67  stp_sequence_t *seq;          /* Sequence (contains the curve data) */
68  double *interval;		/* We allocate an extra slot for the
69				   wrap-around value. */
70
71};
72
73static const char *const stpi_curve_type_names[] =
74  {
75    "linear",
76    "spline",
77  };
78
79static const int stpi_curve_type_count =
80(sizeof(stpi_curve_type_names) / sizeof(const char *));
81
82static const char *const stpi_wrap_mode_names[] =
83  {
84    "nowrap",
85    "wrap"
86  };
87
88static const int stpi_wrap_mode_count =
89(sizeof(stpi_wrap_mode_names) / sizeof(const char *));
90
91/*
92 * Get the total number of points in the base sequence class
93 */
94static size_t
95get_real_point_count(const stp_curve_t *curve)
96{
97  if (curve->piecewise)
98    return stp_sequence_get_size(curve->seq) / 2;
99  else
100    return stp_sequence_get_size(curve->seq);
101}
102
103/*
104 * Get the number of points used by the curve (that are visible to the
105 * user).  This is the real point count, but is decreased by 1 if the
106 * curve wraps around.
107 */
108static size_t
109get_point_count(const stp_curve_t *curve)
110{
111  size_t count = get_real_point_count(curve);
112  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
113    count -= 1;
114
115  return count;
116}
117
118static void
119invalidate_auxiliary_data(stp_curve_t *curve)
120{
121  STP_SAFE_FREE(curve->interval);
122}
123
124static void
125clear_curve_data(stp_curve_t *curve)
126{
127  if (curve->seq)
128    stp_sequence_set_size(curve->seq, 0);
129  curve->recompute_interval = 0;
130  invalidate_auxiliary_data(curve);
131}
132
133static void
134compute_linear_deltas(stp_curve_t *curve)
135{
136  int i;
137  size_t delta_count;
138  size_t seq_point_count;
139  const double *data;
140
141  stp_sequence_get_data(curve->seq, &seq_point_count, &data);
142  if (data == NULL)
143    return;
144
145  delta_count = get_real_point_count(curve);
146
147  if (delta_count <= 1) /* No intervals can be computed */
148    return;
149  delta_count--; /* One less than the real point count.  Note size_t
150		    is unsigned. */
151
152  curve->interval = stp_malloc(sizeof(double) * delta_count);
153  for (i = 0; i < delta_count; i++)
154    {
155      if (curve->piecewise)
156	curve->interval[i] = data[(2 * (i + 1)) + 1] - data[(2 * i) + 1];
157      else
158	curve->interval[i] = data[i + 1] - data[i];
159    }
160}
161
162static void
163compute_spline_deltas_piecewise(stp_curve_t *curve)
164{
165  int i;
166  int k;
167  double *u;
168  double *y2;
169  const double *data = NULL;
170  const stp_curve_point_t *dp;
171  size_t point_count;
172  size_t real_point_count;
173  double sig;
174  double p;
175
176  point_count = get_point_count(curve);
177
178  stp_sequence_get_data(curve->seq, &real_point_count, &data);
179  dp = (const stp_curve_point_t *)data;
180  real_point_count = real_point_count / 2;
181
182  u = stp_malloc(sizeof(double) * real_point_count);
183  y2 = stp_malloc(sizeof(double) * real_point_count);
184
185  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
186    {
187      int reps = 3;
188      int count = reps * real_point_count;
189      double *y2a = stp_malloc(sizeof(double) * count);
190      double *ua = stp_malloc(sizeof(double) * count);
191      y2a[0] = 0.0;
192      ua[0] = 0.0;
193      for (i = 1; i < count - 1; i++)
194	{
195	  int im1 = (i - 1) % point_count;
196	  int ia = i % point_count;
197	  int ip1 = (i + 1) % point_count;
198
199	  sig = (dp[ia].x - dp[im1].x) / (dp[ip1].x - dp[im1].x);
200	  p = sig * y2a[im1] + 2.0;
201	  y2a[i] = (sig - 1.0) / p;
202
203	  ua[i] = ((dp[ip1].y - dp[ia].y) / (dp[ip1].x - dp[ia].x)) -
204	    ((dp[ia].y - dp[im1].y) / (dp[ia].x - dp[im1].x));
205	  ua[i] =
206	    (((6.0 * ua[ia]) / (dp[ip1].x - dp[im1].x)) - (sig * ua[im1])) / p;
207	}
208      y2a[count - 1] = 0.0;
209      for (k = count - 2 ; k >= 0; k--)
210	y2a[k] = y2a[k] * y2a[k + 1] + ua[k];
211      memcpy(u, ua + ((reps / 2) * point_count),
212	     sizeof(double) * real_point_count);
213      memcpy(y2, y2a + ((reps / 2) * point_count),
214	     sizeof(double) * real_point_count);
215      stp_free(y2a);
216      stp_free(ua);
217    }
218  else
219    {
220      int count = real_point_count - 1;
221
222      y2[0] = 0;
223      u[0] = 2 * (dp[1].y - dp[0].y);
224      for (i = 1; i < count; i++)
225	{
226	  int im1 = (i - 1);
227	  int ip1 = (i + 1);
228
229	  sig = (dp[i].x - dp[im1].x) / (dp[ip1].x - dp[im1].x);
230	  p = sig * y2[im1] + 2.0;
231	  y2[i] = (sig - 1.0) / p;
232
233	  u[i] = ((dp[ip1].y - dp[i].y) / (dp[ip1].x - dp[i].x)) -
234	    ((dp[i].y - dp[im1].y) / (dp[i].x - dp[im1].x));
235	  u[i] =
236	    (((6.0 * u[i]) / (dp[ip1].x - dp[im1].x)) - (sig * u[im1])) / p;
237	  stp_deprintf(STP_DBG_CURVE,
238		       "%d sig %f p %f y2 %f u %f x %f %f %f y %f %f %f\n",
239		       i, sig, p, y2[i], u[i],
240		       dp[im1].x, dp[i].x, dp[ip1].x,
241		       dp[im1].y, dp[i].y, dp[ip1].y);
242	}
243      y2[count] = 0.0;
244      u[count] = 0.0;
245      for (k = real_point_count - 2; k >= 0; k--)
246	y2[k] = y2[k] * y2[k + 1] + u[k];
247    }
248
249  curve->interval = y2;
250  stp_free(u);
251}
252
253static void
254compute_spline_deltas_dense(stp_curve_t *curve)
255{
256  int i;
257  int k;
258  double *u;
259  double *y2;
260  const double *y;
261  size_t point_count;
262  size_t real_point_count;
263  double sig;
264  double p;
265
266  point_count = get_point_count(curve);
267
268  stp_sequence_get_data(curve->seq, &real_point_count, &y);
269  u = stp_malloc(sizeof(double) * real_point_count);
270  y2 = stp_malloc(sizeof(double) * real_point_count);
271
272  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
273    {
274      int reps = 3;
275      int count = reps * real_point_count;
276      double *y2a = stp_malloc(sizeof(double) * count);
277      double *ua = stp_malloc(sizeof(double) * count);
278      y2a[0] = 0.0;
279      ua[0] = 0.0;
280      for (i = 1; i < count - 1; i++)
281	{
282	  int im1 = (i - 1);
283	  int ip1 = (i + 1);
284	  int im1a = im1 % point_count;
285	  int ia = i % point_count;
286	  int ip1a = ip1 % point_count;
287
288	  sig = (i - im1) / (ip1 - im1);
289	  p = sig * y2a[im1] + 2.0;
290	  y2a[i] = (sig - 1.0) / p;
291
292	  ua[i] = y[ip1a] - 2 * y[ia] + y[im1a];
293	  ua[i] = 3.0 * ua[i] - sig * ua[im1] / p;
294	}
295      y2a[count - 1] = 0.0;
296      for (k = count - 2 ; k >= 0; k--)
297	y2a[k] = y2a[k] * y2a[k + 1] + ua[k];
298      memcpy(u, ua + ((reps / 2) * point_count),
299	     sizeof(double) * real_point_count);
300      memcpy(y2, y2a + ((reps / 2) * point_count),
301	     sizeof(double) * real_point_count);
302      stp_free(y2a);
303      stp_free(ua);
304    }
305  else
306    {
307      int count = real_point_count - 1;
308
309      y2[0] = 0;
310      u[0] = 2 * (y[1] - y[0]);
311      for (i = 1; i < count; i++)
312	{
313	  int im1 = (i - 1);
314	  int ip1 = (i + 1);
315
316	  sig = (i - im1) / (ip1 - im1);
317	  p = sig * y2[im1] + 2.0;
318	  y2[i] = (sig - 1.0) / p;
319
320	  u[i] = y[ip1] - 2 * y[i] + y[im1];
321	  u[i] = 3.0 * u[i] - sig * u[im1] / p;
322	}
323
324      u[count] = 2 * (y[count] - y[count - 1]);
325      y2[count] = 0.0;
326
327      u[count] = 0.0;
328      for (k = real_point_count - 2; k >= 0; k--)
329	y2[k] = y2[k] * y2[k + 1] + u[k];
330    }
331
332  curve->interval = y2;
333  stp_free(u);
334}
335
336static void
337compute_spline_deltas(stp_curve_t *curve)
338{
339  if (curve->piecewise)
340    compute_spline_deltas_piecewise(curve);
341  else
342    compute_spline_deltas_dense(curve);
343}
344
345/*
346 * Recompute the delta values for interpolation.
347 * When we actually do support spline curves, this routine will
348 * compute the second derivatives for that purpose, too.
349 */
350static void
351compute_intervals(stp_curve_t *curve)
352{
353  if (curve->interval)
354    {
355      stp_free(curve->interval);
356      curve->interval = NULL;
357    }
358  if (stp_sequence_get_size(curve->seq) > 0)
359    {
360      switch (curve->curve_type)
361	{
362	case STP_CURVE_TYPE_SPLINE:
363	  compute_spline_deltas(curve);
364	  break;
365	case STP_CURVE_TYPE_LINEAR:
366	  compute_linear_deltas(curve);
367	  break;
368	}
369    }
370  curve->recompute_interval = 0;
371}
372
373static int
374stpi_curve_set_points(stp_curve_t *curve, size_t points)
375{
376  if (points < 2)
377    return 0;
378  if (points > curve_point_limit ||
379      (curve->wrap_mode == STP_CURVE_WRAP_AROUND &&
380       points > curve_point_limit - 1))
381    return 0;
382  clear_curve_data(curve);
383  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
384    points++;
385  if (curve->piecewise)
386    points *= 2;
387  if ((stp_sequence_set_size(curve->seq, points)) == 0)
388    return 0;
389  return 1;
390}
391
392/*
393 * Create a default curve
394 */
395static void
396stpi_curve_ctor(stp_curve_t *curve, stp_curve_wrap_mode_t wrap_mode)
397{
398  curve->seq = stp_sequence_create();
399  stp_sequence_set_bounds(curve->seq, 0.0, 1.0);
400  curve->curve_type = STP_CURVE_TYPE_LINEAR;
401  curve->wrap_mode = wrap_mode;
402  curve->piecewise = 0;
403  stpi_curve_set_points(curve, 2);
404  curve->recompute_interval = 1;
405  if (wrap_mode == STP_CURVE_WRAP_NONE)
406    curve->gamma = 1.0;
407  stp_sequence_set_point(curve->seq, 0, 0);
408  stp_sequence_set_point(curve->seq, 1, 1);
409}
410
411stp_curve_t *
412stp_curve_create(stp_curve_wrap_mode_t wrap_mode)
413{
414  stp_curve_t *ret;
415  if (wrap_mode != STP_CURVE_WRAP_NONE && wrap_mode != STP_CURVE_WRAP_AROUND)
416    return NULL;
417  ret = stp_zalloc(sizeof(stp_curve_t));
418  stpi_curve_ctor(ret, wrap_mode);
419  return ret;
420}
421
422static void
423curve_dtor(stp_curve_t *curve)
424{
425  CHECK_CURVE(curve);
426  clear_curve_data(curve);
427  if (curve->seq)
428    stp_sequence_destroy(curve->seq);
429  memset(curve, 0, sizeof(stp_curve_t));
430  curve->curve_type = -1;
431}
432
433void
434stp_curve_destroy(stp_curve_t *curve)
435{
436  curve_dtor(curve);
437  stp_free(curve);
438}
439
440void
441stp_curve_copy(stp_curve_t *dest, const stp_curve_t *source)
442{
443  CHECK_CURVE(dest);
444  CHECK_CURVE(source);
445  curve_dtor(dest);
446  dest->curve_type = source->curve_type;
447  dest->wrap_mode = source->wrap_mode;
448  dest->gamma = source->gamma;
449  dest->seq = stp_sequence_create_copy(source->seq);
450  dest->piecewise = source->piecewise;
451  dest->recompute_interval = 1;
452}
453
454stp_curve_t *
455stp_curve_create_copy(const stp_curve_t *curve)
456{
457  stp_curve_t *ret;
458  CHECK_CURVE(curve);
459  ret = stp_curve_create(curve->wrap_mode);
460  stp_curve_copy(ret, curve);
461  return ret;
462}
463
464void
465stp_curve_reverse(stp_curve_t *dest, const stp_curve_t *source)
466{
467  CHECK_CURVE(dest);
468  CHECK_CURVE(source);
469  curve_dtor(dest);
470  dest->curve_type = source->curve_type;
471  dest->wrap_mode = source->wrap_mode;
472  dest->gamma = source->gamma;
473  if (source->piecewise)
474    {
475      const double *source_data;
476      size_t size;
477      double *new_data;
478      int i;
479      stp_sequence_get_data(source->seq, &size, &source_data);
480      new_data = stp_malloc(sizeof(double) * size);
481      for (i = 0; i < size; i += 2)
482	{
483	  int j = size - i - 2;
484	  new_data[i] = 1.0 - source_data[j];
485	  new_data[i + 1] = source_data[j + 1];
486	}
487      dest->seq = stp_sequence_create();
488      stp_sequence_set_data(dest->seq, size, new_data);
489      stp_free(new_data);
490    }
491  else
492      dest->seq = stp_sequence_create_reverse(source->seq);
493  dest->piecewise = source->piecewise;
494  dest->recompute_interval = 1;
495}
496
497stp_curve_t *
498stp_curve_create_reverse(const stp_curve_t *curve)
499{
500  stp_curve_t *ret;
501  CHECK_CURVE(curve);
502  ret = stp_curve_create(curve->wrap_mode);
503  stp_curve_reverse(ret, curve);
504  return ret;
505}
506
507int
508stp_curve_set_bounds(stp_curve_t *curve, double low, double high)
509{
510  CHECK_CURVE(curve);
511  return stp_sequence_set_bounds(curve->seq, low, high);
512}
513
514void
515stp_curve_get_bounds(const stp_curve_t *curve, double *low, double *high)
516{
517  CHECK_CURVE(curve);
518  stp_sequence_get_bounds(curve->seq, low, high);
519}
520
521/*
522 * Find the minimum and maximum points on the curve.  This does not
523 * attempt to find the minimum and maximum interpolations; with cubic
524 * splines these could exceed the boundaries.  That's OK; the interpolation
525 * code will clip them to the bounds.
526 */
527void
528stp_curve_get_range(const stp_curve_t *curve, double *low, double *high)
529{
530  CHECK_CURVE(curve);
531  stp_sequence_get_range(curve->seq, low, high);
532}
533
534size_t
535stp_curve_count_points(const stp_curve_t *curve)
536{
537  CHECK_CURVE(curve);
538  return get_point_count(curve);
539}
540
541stp_curve_wrap_mode_t
542stp_curve_get_wrap(const stp_curve_t *curve)
543{
544  CHECK_CURVE(curve);
545  return curve->wrap_mode;
546}
547
548int
549stp_curve_is_piecewise(const stp_curve_t *curve)
550{
551  CHECK_CURVE(curve);
552  return curve->piecewise;
553}
554
555int
556stp_curve_set_interpolation_type(stp_curve_t *curve, stp_curve_type_t itype)
557{
558  CHECK_CURVE(curve);
559  if (itype < 0 || itype >= stpi_curve_type_count)
560    return 0;
561  curve->curve_type = itype;
562  return 1;
563}
564
565stp_curve_type_t
566stp_curve_get_interpolation_type(const stp_curve_t *curve)
567{
568  CHECK_CURVE(curve);
569  return curve->curve_type;
570}
571
572int
573stp_curve_set_gamma(stp_curve_t *curve, double fgamma)
574{
575  CHECK_CURVE(curve);
576  if (curve->wrap_mode || ! isfinite(fgamma) || fgamma == 0.0)
577    return 0;
578  clear_curve_data(curve);
579  curve->gamma = fgamma;
580  stp_curve_resample(curve, 2);
581  return 1;
582}
583
584double
585stp_curve_get_gamma(const stp_curve_t *curve)
586{
587  CHECK_CURVE(curve);
588  return curve->gamma;
589}
590
591int
592stp_curve_set_data(stp_curve_t *curve, size_t count, const double *data)
593{
594  size_t i;
595  size_t real_count = count;
596  double low, high;
597  CHECK_CURVE(curve);
598  if (count < 2)
599    return 0;
600  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
601    real_count++;
602  if (real_count > curve_point_limit)
603    return 0;
604
605  /* Validate the data before we commit to it. */
606  stp_sequence_get_bounds(curve->seq, &low, &high);
607  for (i = 0; i < count; i++)
608    if (! isfinite(data[i]) || data[i] < low || data[i] > high)
609      {
610	stp_deprintf(STP_DBG_CURVE_ERRORS,
611		     "stp_curve_set_data: datum out of bounds: "
612		     "%g (require %g <= x <= %g), n = %ld\n",
613		     data[i], low, high, (long)i);
614	return 0;
615      }
616  /* Allocate sequence; also accounts for WRAP_MODE */
617  stpi_curve_set_points(curve, count);
618  curve->gamma = 0.0;
619  stp_sequence_set_subrange(curve->seq, 0, count, data);
620
621  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
622    stp_sequence_set_point(curve->seq, count, data[0]);
623  curve->recompute_interval = 1;
624  curve->piecewise = 0;
625
626  return 1;
627}
628
629int
630stp_curve_set_data_points(stp_curve_t *curve, size_t count,
631			  const stp_curve_point_t *data)
632{
633  size_t i;
634  size_t real_count = count;
635  double low, high;
636  double last_x = -1;
637  CHECK_CURVE(curve);
638  if (count < 2)
639    {
640      stp_deprintf(STP_DBG_CURVE_ERRORS,
641		   "stp_curve_set_data_points: too few points %ld\n", (long)count);
642      return 0;
643    }
644  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
645    real_count++;
646  if (real_count > curve_point_limit)
647    {
648      stp_deprintf(STP_DBG_CURVE_ERRORS,
649		   "stp_curve_set_data_points: too many points %ld\n",
650		   (long)real_count);
651      return 0;
652    }
653
654  /* Validate the data before we commit to it. */
655  stp_sequence_get_bounds(curve->seq, &low, &high);
656  for (i = 0; i < count; i++)
657    {
658      if (! isfinite(data[i].y) || data[i].y < low || data[i].y > high)
659	{
660	  stp_deprintf(STP_DBG_CURVE_ERRORS,
661		       "stp_curve_set_data_points: datum out of bounds: "
662		       "%g (require %g <= x <= %g), n = %ld\n",
663		       data[i].y, low, high, (long)i);
664	  return 0;
665	}
666      if (i == 0 && data[i].x != 0.0)
667	{
668	  stp_deprintf(STP_DBG_CURVE_ERRORS,
669		       "stp_curve_set_data_points: first point must have x=0\n");
670	  return 0;
671	}
672      if (curve->wrap_mode == STP_CURVE_WRAP_NONE && i == count - 1 &&
673	  data[i].x != 1.0)
674	{
675	  stp_deprintf(STP_DBG_CURVE_ERRORS,
676		       "stp_curve_set_data_points: last point must have x=1\n");
677	  return 0;
678	}
679      if (curve->wrap_mode == STP_CURVE_WRAP_AROUND &&
680	  data[i].x >= 1.0 - .000001)
681	{
682	  stp_deprintf(STP_DBG_CURVE_ERRORS,
683		       "stp_curve_set_data_points: horizontal value must "
684		       "not exceed .99999\n");
685	  return 0;
686	}
687      if (data[i].x < 0 || data[i].x > 1)
688	{
689	  stp_deprintf(STP_DBG_CURVE_ERRORS,
690		       "stp_curve_set_data_points: horizontal position out of bounds: "
691		       "%g, n = %ld\n",
692		       data[i].x, (long)i);
693	  return 0;
694	}
695      if (data[i].x - .000001 < last_x)
696	{
697	  stp_deprintf(STP_DBG_CURVE_ERRORS,
698		       "stp_curve_set_data_points: horizontal position must "
699		       "exceed previous position by .000001: %g, %g, n = %ld\n",
700		       data[i].x, last_x, (long)i);
701	  return 0;
702	}
703      last_x = data[i].x;
704    }
705  /* Allocate sequence; also accounts for WRAP_MODE */
706  curve->piecewise = 1;
707  stpi_curve_set_points(curve, count);
708  curve->gamma = 0.0;
709  stp_sequence_set_subrange(curve->seq, 0, count * 2, (const double *) data);
710
711  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
712    {
713      stp_sequence_set_point(curve->seq, count * 2, data[0].x);
714      stp_sequence_set_point(curve->seq, count * 2 + 1, data[0].y);
715    }
716  curve->recompute_interval = 1;
717
718  return 1;
719}
720
721
722/*
723 * Note that we return a pointer to the raw data here.
724 * A lot of operations change the data vector, that's why we don't
725 * guarantee it across non-const calls.
726 */
727const double *
728stp_curve_get_data(const stp_curve_t *curve, size_t *count)
729{
730  const double *ret;
731  CHECK_CURVE(curve);
732  if (curve->piecewise)
733    return NULL;
734  stp_sequence_get_data(curve->seq, count, &ret);
735  *count = get_point_count(curve);
736  return ret;
737}
738
739const stp_curve_point_t *
740stp_curve_get_data_points(const stp_curve_t *curve, size_t *count)
741{
742  const double *ret;
743  CHECK_CURVE(curve);
744  if (!curve->piecewise)
745    return NULL;
746  stp_sequence_get_data(curve->seq, count, &ret);
747  *count = get_point_count(curve);
748  return (const stp_curve_point_t *) ret;
749}
750
751static const double *
752stpi_curve_get_data_internal(const stp_curve_t *curve, size_t *count)
753{
754  const double *ret;
755  CHECK_CURVE(curve);
756  stp_sequence_get_data(curve->seq, count, &ret);
757  *count = get_point_count(curve);
758  if (curve->piecewise)
759    *count *= 2;
760  return ret;
761}
762
763
764/* "Overloaded" functions */
765
766#define DEFINE_DATA_SETTER(t, name)                                        \
767int                                                                        \
768stp_curve_set_##name##_data(stp_curve_t *curve, size_t count, const t *data) \
769{                                                                          \
770  double *tmp_data;                                                        \
771  size_t i;                                                                \
772  int status;                                                              \
773  size_t real_count = count;                                               \
774                                                                           \
775  CHECK_CURVE(curve);                                                      \
776  if (count < 2)                                                           \
777    return 0;                                                              \
778  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)                           \
779    real_count++;                                                          \
780  if (real_count > curve_point_limit)                                      \
781    return 0;                                                              \
782  tmp_data = stp_malloc(count * sizeof(double));                           \
783  for (i = 0; i < count; i++)                                              \
784    tmp_data[i] = (double) data[i];                                        \
785  status = stp_curve_set_data(curve, count, tmp_data);                     \
786  stp_free(tmp_data);                                                      \
787  return status;                                                           \
788 }
789
790DEFINE_DATA_SETTER(float, float)
791DEFINE_DATA_SETTER(long, long)
792DEFINE_DATA_SETTER(unsigned long, ulong)
793DEFINE_DATA_SETTER(int, int)
794DEFINE_DATA_SETTER(unsigned int, uint)
795DEFINE_DATA_SETTER(short, short)
796DEFINE_DATA_SETTER(unsigned short, ushort)
797
798
799#define DEFINE_DATA_ACCESSOR(t, name)					\
800const t *								\
801stp_curve_get_##name##_data(const stp_curve_t *curve, size_t *count)    \
802{									\
803  if (curve->piecewise)							\
804    return 0;								\
805  return stp_sequence_get_##name##_data(curve->seq, count);		\
806}
807
808DEFINE_DATA_ACCESSOR(float, float)
809DEFINE_DATA_ACCESSOR(long, long)
810DEFINE_DATA_ACCESSOR(unsigned long, ulong)
811DEFINE_DATA_ACCESSOR(int, int)
812DEFINE_DATA_ACCESSOR(unsigned int, uint)
813DEFINE_DATA_ACCESSOR(short, short)
814DEFINE_DATA_ACCESSOR(unsigned short, ushort)
815
816
817stp_curve_t *
818stp_curve_get_subrange(const stp_curve_t *curve, size_t start, size_t count)
819{
820  stp_curve_t *retval;
821  size_t ncount;
822  double blo, bhi;
823  const double *data;
824  if (start + count > stp_curve_count_points(curve) || count < 2)
825    return NULL;
826  if (curve->piecewise)
827    return NULL;
828  retval = stp_curve_create(STP_CURVE_WRAP_NONE);
829  stp_curve_get_bounds(curve, &blo, &bhi);
830  stp_curve_set_bounds(retval, blo, bhi);
831  data = stp_curve_get_data(curve, &ncount);
832  if (! stp_curve_set_data(retval, count, data + start))
833    {
834      stp_curve_destroy(retval);
835      return NULL;
836    }
837  return retval;
838}
839
840int
841stp_curve_set_subrange(stp_curve_t *curve, const stp_curve_t *range,
842		       size_t start)
843{
844  double blo, bhi;
845  double rlo, rhi;
846  const double *data;
847  size_t count;
848  CHECK_CURVE(curve);
849  if (start + stp_curve_count_points(range) > stp_curve_count_points(curve))
850    return 0;
851  if (curve->piecewise)
852    return 0;
853  stp_sequence_get_bounds(curve->seq, &blo, &bhi);
854  stp_sequence_get_range(curve->seq, &rlo, &rhi);
855  if (rlo < blo || rhi > bhi)
856    return 0;
857  stp_sequence_get_data(range->seq, &count, &data);
858  curve->recompute_interval = 1;
859  curve->gamma = 0.0;
860  invalidate_auxiliary_data(curve);
861  stp_sequence_set_subrange(curve->seq, start, stp_curve_count_points(range),
862			    data);
863  return 1;
864}
865
866
867int
868stp_curve_set_point(stp_curve_t *curve, size_t where, double data)
869{
870  CHECK_CURVE(curve);
871  if (where >= get_point_count(curve))
872    return 0;
873  curve->gamma = 0.0;
874
875  if (curve->piecewise)
876    return 0;
877  if ((stp_sequence_set_point(curve->seq, where, data)) == 0)
878    return 0;
879  if (where == 0 && curve->wrap_mode == STP_CURVE_WRAP_AROUND)
880    if ((stp_sequence_set_point(curve->seq,
881				get_point_count(curve), data)) == 0)
882      return 0;
883  invalidate_auxiliary_data(curve);
884  return 1;
885}
886
887int
888stp_curve_get_point(const stp_curve_t *curve, size_t where, double *data)
889{
890  CHECK_CURVE(curve);
891  if (where >= get_point_count(curve))
892    return 0;
893  if (curve->piecewise)
894    return 0;
895  return stp_sequence_get_point(curve->seq, where, data);
896}
897
898const stp_sequence_t *
899stp_curve_get_sequence(const stp_curve_t *curve)
900{
901  CHECK_CURVE(curve);
902  if (curve->piecewise)
903    return NULL;
904  return curve->seq;
905}
906
907int
908stp_curve_rescale(stp_curve_t *curve, double scale,
909		  stp_curve_compose_t mode, stp_curve_bounds_t bounds_mode)
910{
911  size_t real_point_count;
912  int i;
913  double nblo;
914  double nbhi;
915  size_t count;
916
917  CHECK_CURVE(curve);
918
919  real_point_count = get_real_point_count(curve);
920
921  stp_sequence_get_bounds(curve->seq, &nblo, &nbhi);
922  if (bounds_mode == STP_CURVE_BOUNDS_RESCALE)
923    {
924      switch (mode)
925	{
926	case STP_CURVE_COMPOSE_ADD:
927	  nblo += scale;
928	  nbhi += scale;
929	  break;
930	case STP_CURVE_COMPOSE_MULTIPLY:
931	  if (scale < 0)
932	    {
933	      double tmp = nblo * scale;
934	      nblo = nbhi * scale;
935	      nbhi = tmp;
936	    }
937	  else
938	    {
939	      nblo *= scale;
940	      nbhi *= scale;
941	    }
942	  break;
943	case STP_CURVE_COMPOSE_EXPONENTIATE:
944	  if (scale == 0.0)
945	    return 0;
946	  if (nblo < 0)
947	    return 0;
948	  nblo = pow(nblo, scale);
949	  nbhi = pow(nbhi, scale);
950	  break;
951	default:
952	  return 0;
953	}
954    }
955
956  if (! isfinite(nbhi) || ! isfinite(nblo))
957    return 0;
958
959  count = get_point_count(curve);
960  if (count)
961    {
962      double *tmp;
963      size_t scount;
964      int stride = 1;
965      int offset = 0;
966      const double *data;
967      if (curve->piecewise)
968	{
969	  stride = 2;
970	  offset = 1;
971	}
972      stp_sequence_get_data(curve->seq, &scount, &data);
973      tmp = stp_malloc(sizeof(double) * scount);
974      memcpy(tmp, data, scount * sizeof(double));
975      for (i = offset; i < scount; i += stride)
976	{
977	  switch (mode)
978	    {
979	    case STP_CURVE_COMPOSE_ADD:
980	      tmp[i] = tmp[i] + scale;
981	      break;
982	    case STP_CURVE_COMPOSE_MULTIPLY:
983	      tmp[i] = tmp[i] * scale;
984	      break;
985	    case STP_CURVE_COMPOSE_EXPONENTIATE:
986	      tmp[i] = pow(tmp[i], scale);
987	      break;
988	    }
989	  if (tmp[i] > nbhi || tmp[i] < nblo)
990	    {
991	      if (bounds_mode == STP_CURVE_BOUNDS_ERROR)
992		{
993		  stp_free(tmp);
994		  return(0);
995		}
996	      else if (tmp[i] > nbhi)
997		tmp[i] = nbhi;
998	      else
999		tmp[i] = nblo;
1000	    }
1001	}
1002      stp_sequence_set_bounds(curve->seq, nblo, nbhi);
1003      curve->gamma = 0.0;
1004      stpi_curve_set_points(curve, count);
1005      stp_sequence_set_subrange(curve->seq, 0, scount, tmp);
1006      stp_free(tmp);
1007      curve->recompute_interval = 1;
1008      invalidate_auxiliary_data(curve);
1009    }
1010  return 1;
1011}
1012
1013static int
1014stpi_curve_check_parameters(stp_curve_t *curve, size_t points)
1015{
1016  double blo, bhi;
1017  if (curve->gamma && curve->wrap_mode)
1018    {
1019      stp_deprintf(STP_DBG_CURVE_ERRORS,
1020		   "curve sets both gamma and wrap_mode\n");
1021      return 0;
1022    }
1023  stp_sequence_get_bounds(curve->seq, &blo, &bhi);
1024  if (blo > bhi)
1025    {
1026      stp_deprintf(STP_DBG_CURVE_ERRORS,
1027		   "curve low bound is greater than high bound\n");
1028      return 0;
1029    }
1030  return 1;
1031}
1032
1033static inline double
1034interpolate_gamma_internal(const stp_curve_t *curve, double where)
1035{
1036  double fgamma = curve->gamma;
1037  double blo, bhi;
1038  size_t real_point_count;
1039
1040  real_point_count = get_real_point_count(curve);;
1041
1042  if (real_point_count)
1043    where /= (real_point_count - 1);
1044  if (fgamma < 0)
1045    {
1046      where = 1.0 - where;
1047      fgamma = -fgamma;
1048    }
1049  stp_sequence_get_bounds(curve->seq, &blo, &bhi);
1050  stp_deprintf(STP_DBG_CURVE,
1051	       "interpolate_gamma %f %f %f %f %f\n", where, fgamma,
1052	       blo, bhi, pow(where, fgamma));
1053  return blo + (bhi - blo) * pow(where, fgamma);
1054}
1055
1056static inline double
1057do_interpolate_spline(double low, double high, double frac,
1058		      double interval_low, double interval_high,
1059		      double x_interval)
1060{
1061  double a = 1.0 - frac;
1062  double b = frac;
1063  double retval =
1064    ((a * a * a - a) * interval_low) + ((b * b * b - b) * interval_high);
1065  retval = retval * x_interval * x_interval / 6;
1066  retval += (a * low) + (b * high);
1067  return retval;
1068}
1069
1070static inline double
1071interpolate_point_internal(stp_curve_t *curve, double where)
1072{
1073  int integer = where;
1074  double frac = where - (double) integer;
1075  double bhi, blo;
1076
1077  if (frac == 0.0)
1078    {
1079      double val;
1080      if ((stp_sequence_get_point(curve->seq, integer, &val)) == 0)
1081	return HUGE_VAL; /* Infinity */
1082      return val;
1083    }
1084  if (curve->recompute_interval)
1085    compute_intervals(curve);
1086  if (curve->curve_type == STP_CURVE_TYPE_LINEAR)
1087    {
1088      double val;
1089      if ((stp_sequence_get_point(curve->seq, integer, &val)) == 0)
1090	return HUGE_VAL; /* Infinity */
1091      return val + frac * curve->interval[integer];
1092    }
1093  else
1094    {
1095      size_t point_count;
1096      double ival, ip1val;
1097      double retval;
1098      int i = integer;
1099      int ip1 = integer + 1;
1100
1101      point_count = get_point_count(curve);
1102
1103      if (ip1 >= point_count)
1104	ip1 -= point_count;
1105
1106      if ((stp_sequence_get_point(curve->seq, i, &ival)) == 0 ||
1107	  (stp_sequence_get_point(curve->seq, ip1, &ip1val)) == 0)
1108	return HUGE_VAL; /* Infinity */
1109
1110      retval = do_interpolate_spline(ival, ip1val, frac, curve->interval[i],
1111				     curve->interval[ip1], 1.0);
1112
1113      stp_sequence_get_bounds(curve->seq, &blo, &bhi);
1114      if (retval > bhi)
1115	retval = bhi;
1116      if (retval < blo)
1117	retval = blo;
1118      return retval;
1119    }
1120}
1121
1122int
1123stp_curve_interpolate_value(const stp_curve_t *curve, double where,
1124			    double *result)
1125{
1126  size_t limit;
1127
1128  CHECK_CURVE(curve);
1129  if (curve->piecewise)
1130    return 0;
1131
1132  limit = get_real_point_count(curve);
1133
1134  if (where < 0 || where > limit)
1135    return 0;
1136  if (curve->gamma)	/* this means a pure gamma curve */
1137    *result = interpolate_gamma_internal(curve, where);
1138  else
1139    *result = interpolate_point_internal((stp_curve_t *) curve, where);
1140  return 1;
1141}
1142
1143int
1144stp_curve_resample(stp_curve_t *curve, size_t points)
1145{
1146  size_t limit = points;
1147  size_t old;
1148  size_t i;
1149  double *new_vec;
1150
1151  CHECK_CURVE(curve);
1152
1153  if (points == get_point_count(curve) && curve->seq && !(curve->piecewise))
1154    return 1;
1155
1156  if (points < 2)
1157    return 1;
1158
1159  if (curve->wrap_mode == STP_CURVE_WRAP_AROUND)
1160    limit++;
1161  if (limit > curve_point_limit)
1162    return 0;
1163  old = get_real_point_count(curve);
1164  if (old)
1165    old--;
1166  if (!old)
1167    old = 1;
1168
1169  new_vec = stp_malloc(sizeof(double) * limit);
1170
1171  /*
1172   * Be very careful how we calculate the location along the scale!
1173   * If we're not careful how we do it, we might get a small roundoff
1174   * error
1175   */
1176  if (curve->piecewise)
1177    {
1178      double blo, bhi;
1179      int curpos = 0;
1180      stp_sequence_get_bounds(curve->seq, &blo, &bhi);
1181      if (curve->recompute_interval)
1182	compute_intervals(curve);
1183      for (i = 0; i < old; i++)
1184	{
1185	  double low;
1186	  double high;
1187	  double low_y;
1188	  double high_y;
1189	  double x_delta;
1190	  if (!stp_sequence_get_point(curve->seq, i * 2, &low))
1191	    {
1192	      stp_free(new_vec);
1193	      return 0;
1194	    }
1195	  if (i == old - 1)
1196	    high = 1.0;
1197	  else if (!stp_sequence_get_point(curve->seq, ((i + 1) * 2), &high))
1198	    {
1199	      stp_free(new_vec);
1200	      return 0;
1201	    }
1202	  if (!stp_sequence_get_point(curve->seq, (i * 2) + 1, &low_y))
1203	    {
1204	      stp_free(new_vec);
1205	      return 0;
1206	    }
1207	  if (!stp_sequence_get_point(curve->seq, ((i + 1) * 2) + 1, &high_y))
1208	    {
1209	      stp_free(new_vec);
1210	      return 0;
1211	    }
1212	  stp_deprintf(STP_DBG_CURVE,
1213		       "Filling slots at %ld %d: %f %f  %f %f  %ld\n",
1214		       (long)i,curpos, high, low, high_y, low_y, (long)limit);
1215	  x_delta = high - low;
1216	  high *= (limit - 1);
1217	  low *= (limit - 1);
1218	  while (curpos <= high)
1219	    {
1220	      double frac = (curpos - low) / (high - low);
1221	      if (curve->curve_type == STP_CURVE_TYPE_LINEAR)
1222		new_vec[curpos] = low_y + frac * (high_y - low_y);
1223	      else
1224		new_vec[curpos] =
1225		  do_interpolate_spline(low_y, high_y, frac,
1226					curve->interval[i],
1227					curve->interval[i + 1],
1228					x_delta);
1229	      if (new_vec[curpos] < blo)
1230		new_vec[curpos] = blo;
1231	      if (new_vec[curpos] > bhi)
1232		new_vec[curpos] = bhi;
1233	      stp_deprintf(STP_DBG_CURVE,
1234			   "  Filling slot %d %f %f\n",
1235			   curpos, frac, new_vec[curpos]);
1236	      curpos++;
1237	    }
1238	}
1239      curve->piecewise = 0;
1240    }
1241  else
1242    {
1243      for (i = 0; i < limit; i++)
1244	if (curve->gamma)
1245	  new_vec[i] =
1246	    interpolate_gamma_internal(curve, ((double) i * (double) old /
1247					       (double) (limit - 1)));
1248	else
1249	  new_vec[i] =
1250	    interpolate_point_internal(curve, ((double) i * (double) old /
1251					       (double) (limit - 1)));
1252    }
1253  stpi_curve_set_points(curve, points);
1254  stp_sequence_set_subrange(curve->seq, 0, limit, new_vec);
1255  curve->recompute_interval = 1;
1256  stp_free(new_vec);
1257  return 1;
1258}
1259
1260static unsigned
1261gcd(unsigned a, unsigned b)
1262{
1263  unsigned tmp;
1264  if (b > a)
1265    {
1266      tmp = a;
1267      a = b;
1268      b = tmp;
1269    }
1270  while (1)
1271    {
1272      tmp = a % b;
1273      if (tmp == 0)
1274	return b;
1275      a = b;
1276      b = tmp;
1277    }
1278}
1279
1280static unsigned
1281lcm(unsigned a, unsigned b)
1282{
1283  if (a == b)
1284    return a;
1285  else if (a * b == 0)
1286    return a > b ? a : b;
1287  else
1288    {
1289      double rval = (double) a / gcd(a, b) * b;
1290      if (rval > curve_point_limit)
1291	return curve_point_limit;
1292      else
1293	return rval;
1294    }
1295}
1296
1297static int
1298create_gamma_curve(stp_curve_t **retval, double lo, double hi, double fgamma,
1299		   int points)
1300{
1301  *retval = stp_curve_create(STP_CURVE_WRAP_NONE);
1302  if (stp_curve_set_bounds(*retval, lo, hi) &&
1303      stp_curve_set_gamma(*retval, fgamma) &&
1304      stp_curve_resample(*retval, points))
1305    return 1;
1306  stp_curve_destroy(*retval);
1307  *retval = 0;
1308  return 0;
1309}
1310
1311static int
1312interpolate_points(stp_curve_t *a, stp_curve_t *b,
1313		   stp_curve_compose_t mode,
1314		   int points, double *tmp_data)
1315{
1316  double pa, pb;
1317  int i;
1318  size_t points_a = stp_curve_count_points(a);
1319  size_t points_b = stp_curve_count_points(b);
1320  for (i = 0; i < points; i++)
1321    {
1322      if (!stp_curve_interpolate_value
1323	  (a, (double) i * (points_a - 1) / (points - 1), &pa))
1324	{
1325	  stp_deprintf(STP_DBG_CURVE_ERRORS,
1326		       "interpolate_points: interpolate curve a value failed\n");
1327	  return 0;
1328	}
1329      if (!stp_curve_interpolate_value
1330	  (b, (double) i * (points_b - 1) / (points - 1), &pb))
1331	{
1332	  stp_deprintf(STP_DBG_CURVE_ERRORS,
1333		       "interpolate_points: interpolate curve b value failed\n");
1334	  return 0;
1335	}
1336      if (mode == STP_CURVE_COMPOSE_ADD)
1337	pa += pb;
1338      else
1339	pa *= pb;
1340      if (! isfinite(pa))
1341	{
1342	  stp_deprintf(STP_DBG_CURVE_ERRORS,
1343		       "interpolate_points: interpolated point %lu is invalid\n",
1344		       (unsigned long) i);
1345	  return 0;
1346	}
1347      tmp_data[i] = pa;
1348    }
1349  return 1;
1350}
1351
1352int
1353stp_curve_compose(stp_curve_t **retval,
1354		  stp_curve_t *a, stp_curve_t *b,
1355		  stp_curve_compose_t mode, int points)
1356{
1357  stp_curve_t *ret;
1358  double *tmp_data;
1359  double gamma_a = stp_curve_get_gamma(a);
1360  double gamma_b = stp_curve_get_gamma(b);
1361  unsigned points_a = stp_curve_count_points(a);
1362  unsigned points_b = stp_curve_count_points(b);
1363  double alo, ahi, blo, bhi;
1364
1365  if (a->piecewise && b->piecewise)
1366    return 0;
1367  if (a->piecewise)
1368    {
1369      stp_curve_t *a_save = a;
1370      a = stp_curve_create_copy(a_save);
1371      stp_curve_resample(a, stp_curve_count_points(b));
1372    }
1373  if (b->piecewise)
1374    {
1375      stp_curve_t *b_save = b;
1376      b = stp_curve_create_copy(b_save);
1377      stp_curve_resample(b, stp_curve_count_points(a));
1378    }
1379
1380  if (mode != STP_CURVE_COMPOSE_ADD && mode != STP_CURVE_COMPOSE_MULTIPLY)
1381    return 0;
1382  if (stp_curve_get_wrap(a) != stp_curve_get_wrap(b))
1383    return 0;
1384  stp_curve_get_bounds(a, &alo, &ahi);
1385  stp_curve_get_bounds(b, &blo, &bhi);
1386  if (mode == STP_CURVE_COMPOSE_MULTIPLY && (alo < 0 || blo < 0))
1387    return 0;
1388
1389  if (stp_curve_get_wrap(a) == STP_CURVE_WRAP_AROUND)
1390    {
1391      points_a++;
1392      points_b++;
1393    }
1394  if (points == -1)
1395    {
1396      points = lcm(points_a, points_b);
1397      if (stp_curve_get_wrap(a) == STP_CURVE_WRAP_AROUND)
1398	points--;
1399    }
1400  if (points < 2 || points > curve_point_limit ||
1401      ((stp_curve_get_wrap(a) == STP_CURVE_WRAP_AROUND) &&
1402       points > curve_point_limit - 1))
1403    return 0;
1404
1405  if (gamma_a && gamma_b && gamma_a * gamma_b > 0 &&
1406      mode == STP_CURVE_COMPOSE_MULTIPLY)
1407    return create_gamma_curve(retval, alo * blo, ahi * bhi, gamma_a + gamma_b,
1408			      points);
1409  tmp_data = stp_malloc(sizeof(double) * points);
1410  if (!interpolate_points(a, b, mode, points, tmp_data))
1411    {
1412      stp_free(tmp_data);
1413      return 0;
1414    }
1415  ret = stp_curve_create(stp_curve_get_wrap(a));
1416  if (mode == STP_CURVE_COMPOSE_ADD)
1417    {
1418      stp_curve_rescale(ret, (ahi - alo) + (bhi - blo),
1419			STP_CURVE_COMPOSE_MULTIPLY, STP_CURVE_BOUNDS_RESCALE);
1420      stp_curve_rescale(ret, alo + blo,
1421			STP_CURVE_COMPOSE_ADD, STP_CURVE_BOUNDS_RESCALE);
1422    }
1423  else
1424    {
1425      stp_curve_rescale(ret, (ahi - alo) * (bhi - blo),
1426			STP_CURVE_COMPOSE_MULTIPLY, STP_CURVE_BOUNDS_RESCALE);
1427      stp_curve_rescale(ret, alo * blo,
1428			STP_CURVE_COMPOSE_ADD, STP_CURVE_BOUNDS_RESCALE);
1429    }
1430  if (! stp_curve_set_data(ret, points, tmp_data))
1431    goto bad1;
1432  *retval = ret;
1433  stp_free(tmp_data);
1434  return 1;
1435 bad1:
1436  stp_curve_destroy(ret);
1437  stp_free(tmp_data);
1438  return 0;
1439}
1440
1441
1442stp_curve_t *
1443stp_curve_create_from_xmltree(stp_mxml_node_t *curve)  /* The curve node */
1444{
1445  const char *stmp;                       /* Temporary string */
1446  stp_mxml_node_t *child;                 /* Child sequence node */
1447  stp_curve_t *ret = NULL;                /* Curve to return */
1448  stp_curve_type_t curve_type;            /* Type of curve */
1449  stp_curve_wrap_mode_t wrap_mode;        /* Curve wrap mode */
1450  double fgamma;                          /* Gamma value */
1451  stp_sequence_t *seq = NULL;             /* Sequence data */
1452  double low, high;                       /* Sequence bounds */
1453  int piecewise = 0;
1454
1455  stp_xml_init();
1456  /* Get curve type */
1457  stmp = stp_mxmlElementGetAttr(curve, "type");
1458  if (stmp)
1459    {
1460      if (!strcmp(stmp, "linear"))
1461	  curve_type = STP_CURVE_TYPE_LINEAR;
1462      else if (!strcmp(stmp, "spline"))
1463	  curve_type = STP_CURVE_TYPE_SPLINE;
1464      else
1465	{
1466	  stp_deprintf(STP_DBG_CURVE_ERRORS,
1467		       "stp_curve_create_from_xmltree: %s: \"type\" invalid\n", stmp);
1468	  goto error;
1469	}
1470    }
1471  else
1472    {
1473      stp_deprintf(STP_DBG_CURVE_ERRORS,
1474		   "stp_curve_create_from_xmltree: \"type\" missing\n");
1475      goto error;
1476    }
1477  /* Get curve wrap mode */
1478  stmp = stp_mxmlElementGetAttr(curve, "wrap");
1479  if (stmp)
1480    {
1481      if (!strcmp(stmp, "nowrap"))
1482	wrap_mode = STP_CURVE_WRAP_NONE;
1483      else if (!strcmp(stmp, "wrap"))
1484	{
1485	  wrap_mode = STP_CURVE_WRAP_AROUND;
1486	}
1487      else
1488	{
1489	  stp_deprintf(STP_DBG_CURVE_ERRORS,
1490		       "stp_curve_create_from_xmltree: %s: \"wrap\" invalid\n", stmp);
1491	  goto error;
1492	}
1493    }
1494  else
1495    {
1496      stp_deprintf(STP_DBG_CURVE_ERRORS,
1497		   "stp_curve_create_from_xmltree: \"wrap\" missing\n");
1498      goto error;
1499    }
1500  /* Get curve gamma */
1501  stmp = stp_mxmlElementGetAttr(curve, "gamma");
1502  if (stmp)
1503    {
1504      fgamma = stp_xmlstrtod(stmp);
1505    }
1506  else
1507    {
1508      stp_deprintf(STP_DBG_CURVE_ERRORS,
1509		   "stp_curve_create_from_xmltree: \"gamma\" missing\n");
1510      goto error;
1511    }
1512  /* If gamma is set, wrap_mode must be STP_CURVE_WRAP_NONE */
1513  if (fgamma && wrap_mode != STP_CURVE_WRAP_NONE)
1514    {
1515      stp_deprintf(STP_DBG_CURVE_ERRORS, "stp_curve_create_from_xmltree: "
1516		   "gamma set and \"wrap\" is not STP_CURVE_WRAP_NONE\n");
1517      goto error;
1518    }
1519  stmp = stp_mxmlElementGetAttr(curve, "piecewise");
1520  if (stmp && strcmp(stmp, "true") == 0)
1521    piecewise = 1;
1522
1523  /* Set up the curve */
1524  ret = stp_curve_create(wrap_mode);
1525  stp_curve_set_interpolation_type(ret, curve_type);
1526
1527  child = stp_mxmlFindElement(curve, curve, "sequence", NULL, NULL, STP_MXML_DESCEND);
1528  if (child)
1529    seq = stp_sequence_create_from_xmltree(child);
1530
1531  if (seq == NULL)
1532    {
1533      stp_deprintf(STP_DBG_CURVE_ERRORS,
1534		   "stp_curve_create_from_xmltree: sequence read failed\n");
1535      goto error;
1536    }
1537
1538  /* Set curve bounds */
1539  stp_sequence_get_bounds(seq, &low, &high);
1540  stp_curve_set_bounds(ret, low, high);
1541
1542  if (fgamma)
1543    stp_curve_set_gamma(ret, fgamma);
1544  else /* Not a gamma curve, so set points */
1545    {
1546      size_t seq_count;
1547      const double* data;
1548
1549      stp_sequence_get_data(seq, &seq_count, &data);
1550      if (piecewise)
1551	{
1552	  if ((seq_count % 2) != 0)
1553	    {
1554	      stp_deprintf(STP_DBG_CURVE_ERRORS,
1555			   "stp_curve_create_from_xmltree: invalid data count %ld\n",
1556			   (long)seq_count);
1557	      goto error;
1558	    }
1559	  if (stp_curve_set_data_points(ret, seq_count / 2,
1560					(const stp_curve_point_t *) data) == 0)
1561	    {
1562	      stp_deprintf(STP_DBG_CURVE_ERRORS,
1563			   "stp_curve_create_from_xmltree: failed to set curve data points\n");
1564	      goto error;
1565	    }
1566	}
1567      else
1568	{
1569	  if (stp_curve_set_data(ret, seq_count, data) == 0)
1570	    {
1571	      stp_deprintf(STP_DBG_CURVE_ERRORS,
1572			   "stp_curve_create_from_xmltree: failed to set curve data\n");
1573	      goto error;
1574	    }
1575	}
1576    }
1577
1578  if (seq)
1579    {
1580      stp_sequence_destroy(seq);
1581      seq = NULL;
1582    }
1583
1584    /* Validate curve */
1585  if (stpi_curve_check_parameters(ret, stp_curve_count_points(ret)) == 0)
1586    {
1587      stp_deprintf(STP_DBG_CURVE_ERRORS,
1588		   "stp_curve_create_from_xmltree: parameter check failed\n");
1589      goto error;
1590    }
1591
1592  stp_xml_exit();
1593
1594  return ret;
1595
1596 error:
1597  stp_deprintf(STP_DBG_CURVE_ERRORS,
1598	       "stp_curve_create_from_xmltree: error during curve read\n");
1599  if (ret)
1600    stp_curve_destroy(ret);
1601  stp_xml_exit();
1602  return NULL;
1603}
1604
1605
1606stp_mxml_node_t *
1607stp_xmltree_create_from_curve(const stp_curve_t *curve)  /* The curve */
1608{
1609  stp_curve_wrap_mode_t wrapmode;
1610  stp_curve_type_t interptype;
1611  double gammaval, low, high;
1612  stp_sequence_t *seq;
1613
1614  char *cgamma;
1615
1616  stp_mxml_node_t *curvenode = NULL;
1617  stp_mxml_node_t *child = NULL;
1618
1619  stp_xml_init();
1620
1621  /* Get curve details */
1622  wrapmode = stp_curve_get_wrap(curve);
1623  interptype = stp_curve_get_interpolation_type(curve);
1624  gammaval = stp_curve_get_gamma(curve);
1625
1626  if (gammaval && wrapmode != STP_CURVE_WRAP_NONE)
1627    {
1628      stp_deprintf(STP_DBG_CURVE_ERRORS, "stp_xmltree_create_from_curve: "
1629		   "curve sets gamma and wrap_mode is not STP_CURVE_WRAP_NONE\n");
1630      goto error;
1631    }
1632
1633  /* Construct the allocated strings required */
1634  stp_asprintf(&cgamma, "%g", gammaval);
1635
1636  curvenode = stp_mxmlNewElement(NULL, "curve");
1637  stp_mxmlElementSetAttr(curvenode, "wrap", stpi_wrap_mode_names[wrapmode]);
1638  stp_mxmlElementSetAttr(curvenode, "type", stpi_curve_type_names[interptype]);
1639  stp_mxmlElementSetAttr(curvenode, "gamma", cgamma);
1640  if (curve->piecewise)
1641    stp_mxmlElementSetAttr(curvenode, "piecewise", "true");
1642  else
1643    stp_mxmlElementSetAttr(curvenode, "piecewise", "false");
1644
1645  stp_free(cgamma);
1646
1647  seq = stp_sequence_create();
1648  stp_curve_get_bounds(curve, &low, &high);
1649  stp_sequence_set_bounds(seq, low, high);
1650  if (gammaval != 0) /* A gamma curve does not require sequence data */
1651    {
1652      stp_sequence_set_size(seq, 0);
1653    }
1654  else
1655    {
1656      const double *data;
1657      size_t count;
1658      data = stpi_curve_get_data_internal(curve, &count);
1659      stp_sequence_set_data(seq, count, data);
1660    }
1661
1662  child = stp_xmltree_create_from_sequence(seq);
1663
1664  if (seq)
1665    {
1666      stp_sequence_destroy(seq);
1667      seq = NULL;
1668    }
1669
1670  if (child == NULL)
1671    {
1672      stp_deprintf(STP_DBG_CURVE_ERRORS,
1673		   "stp_xmltree_create_from_curve: sequence node is NULL\n");
1674      goto error;
1675    }
1676  stp_mxmlAdd(curvenode, STP_MXML_ADD_AFTER, NULL, child);
1677
1678  stp_xml_exit();
1679
1680  return curvenode;
1681
1682 error:
1683  stp_deprintf(STP_DBG_CURVE_ERRORS,
1684	       "stp_xmltree_create_from_curve: error during xmltree creation\n");
1685  if (curvenode)
1686    stp_mxmlDelete(curvenode);
1687  if (child)
1688    stp_mxmlDelete(child);
1689  stp_xml_exit();
1690
1691  return NULL;
1692}
1693
1694static stp_mxml_node_t *
1695xmldoc_create_from_curve(const stp_curve_t *curve)
1696{
1697  stp_mxml_node_t *xmldoc;
1698  stp_mxml_node_t *rootnode;
1699  stp_mxml_node_t *curvenode;
1700
1701  /* Get curve details */
1702  curvenode = stp_xmltree_create_from_curve(curve);
1703  if (curvenode == NULL)
1704    {
1705      stp_deprintf(STP_DBG_CURVE_ERRORS,
1706		   "xmldoc_create_from_curve: error creating curve node\n");
1707      return NULL;
1708    }
1709  /* Create the XML tree */
1710  xmldoc = stp_xmldoc_create_generic();
1711  if (xmldoc == NULL)
1712    {
1713      stp_deprintf(STP_DBG_CURVE_ERRORS,
1714		   "xmldoc_create_from_curve: error creating XML document\n");
1715      return NULL;
1716    }
1717  rootnode = xmldoc->child;
1718  if (rootnode == NULL)
1719    {
1720      stp_mxmlDelete(xmldoc);
1721      stp_deprintf(STP_DBG_CURVE_ERRORS,
1722		   "xmldoc_create_from_curve: error getting XML document root node\n");
1723      return NULL;
1724    }
1725
1726  stp_mxmlAdd(rootnode, STP_MXML_ADD_AFTER, NULL, curvenode);
1727
1728  return xmldoc;
1729}
1730
1731static int
1732curve_whitespace_callback(stp_mxml_node_t *node, int where)
1733{
1734  if (node->type != STP_MXML_ELEMENT)
1735    return 0;
1736  if (strcasecmp(node->value.element.name, "gutenprint") == 0)
1737    {
1738      switch (where)
1739	{
1740	case STP_MXML_WS_AFTER_OPEN:
1741	case STP_MXML_WS_BEFORE_CLOSE:
1742	case STP_MXML_WS_AFTER_CLOSE:
1743	  return '\n';
1744	case STP_MXML_WS_BEFORE_OPEN:
1745	default:
1746	  return 0;
1747	}
1748    }
1749  else if (strcasecmp(node->value.element.name, "curve") == 0)
1750    {
1751      switch (where)
1752	{
1753	case STP_MXML_WS_AFTER_OPEN:
1754	  return '\n';
1755	case STP_MXML_WS_BEFORE_CLOSE:
1756	case STP_MXML_WS_AFTER_CLOSE:
1757	case STP_MXML_WS_BEFORE_OPEN:
1758	default:
1759	  return 0;
1760	}
1761    }
1762  else if (strcasecmp(node->value.element.name, "sequence") == 0)
1763    {
1764      const char *count;
1765      switch (where)
1766	{
1767	case STP_MXML_WS_BEFORE_CLOSE:
1768	  count = stp_mxmlElementGetAttr(node, "count");
1769	  if (strcmp(count, "0") == 0)
1770	    return 0;
1771	  else
1772	    return '\n';
1773	case STP_MXML_WS_AFTER_OPEN:
1774	case STP_MXML_WS_AFTER_CLOSE:
1775	  return '\n';
1776	case STP_MXML_WS_BEFORE_OPEN:
1777	default:
1778	  return 0;
1779	}
1780    }
1781  else
1782    return 0;
1783}
1784
1785
1786int
1787stp_curve_write(FILE *file, const stp_curve_t *curve)  /* The curve */
1788{
1789  stp_mxml_node_t *xmldoc = NULL;
1790
1791  stp_xml_init();
1792
1793  xmldoc = xmldoc_create_from_curve(curve);
1794  if (xmldoc == NULL)
1795    {
1796      stp_xml_exit();
1797      return 1;
1798    }
1799
1800  stp_mxmlSaveFile(xmldoc, file, curve_whitespace_callback);
1801
1802  if (xmldoc)
1803    stp_mxmlDelete(xmldoc);
1804
1805  stp_xml_exit();
1806
1807  return 0;
1808}
1809
1810char *
1811stp_curve_write_string(const stp_curve_t *curve)  /* The curve */
1812{
1813  stp_mxml_node_t *xmldoc = NULL;
1814  char *retval;
1815
1816  stp_xml_init();
1817
1818  xmldoc = xmldoc_create_from_curve(curve);
1819  if (xmldoc == NULL)
1820    {
1821      stp_xml_exit();
1822      return NULL;
1823    }
1824
1825  retval = stp_mxmlSaveAllocString(xmldoc, curve_whitespace_callback);
1826
1827  if (xmldoc)
1828    stp_mxmlDelete(xmldoc);
1829
1830  stp_xml_exit();
1831
1832  return retval;
1833}
1834
1835static stp_curve_t *
1836xml_doc_get_curve(stp_mxml_node_t *doc)
1837{
1838  stp_mxml_node_t *cur;
1839  stp_mxml_node_t *xmlcurve;
1840  stp_curve_t *curve = NULL;
1841
1842  if (doc == NULL )
1843    {
1844      stp_deprintf(STP_DBG_CURVE_ERRORS,
1845		   "xml_doc_get_curve: XML file not parsed successfully.\n");
1846      return NULL;
1847    }
1848
1849  cur = doc->child;
1850
1851  if (cur == NULL)
1852    {
1853      stp_deprintf(STP_DBG_CURVE_ERRORS,
1854		   "xml_doc_get_curve: empty document\n");
1855      return NULL;
1856    }
1857
1858  xmlcurve = stp_xml_get_node(cur, "gutenprint", "curve", NULL);
1859
1860  if (xmlcurve)
1861    curve = stp_curve_create_from_xmltree(xmlcurve);
1862
1863  return curve;
1864}
1865
1866stp_curve_t *
1867stp_curve_create_from_file(const char* file)
1868{
1869  stp_curve_t *curve = NULL;
1870  stp_mxml_node_t *doc;
1871  FILE *fp = fopen(file, "r");
1872  if (!fp)
1873    {
1874      stp_deprintf(STP_DBG_CURVE_ERRORS,
1875		   "stp_curve_create_from_file: unable to open %s: %s\n",
1876		    file, strerror(errno));
1877      return NULL;
1878    }
1879  stp_deprintf(STP_DBG_XML, "stp_curve_create_from_file: reading `%s'...\n",
1880	       file);
1881
1882  stp_xml_init();
1883
1884  doc = stp_mxmlLoadFile(NULL, fp, STP_MXML_NO_CALLBACK);
1885
1886  curve = xml_doc_get_curve(doc);
1887
1888  if (doc)
1889    stp_mxmlDelete(doc);
1890
1891  stp_xml_exit();
1892  (void) fclose(fp);
1893  return curve;
1894
1895}
1896
1897stp_curve_t *
1898stp_curve_create_from_stream(FILE* fp)
1899{
1900  stp_curve_t *curve = NULL;
1901  stp_mxml_node_t *doc;
1902  stp_deprintf(STP_DBG_XML, "stp_curve_create_from_fp: reading...\n");
1903
1904  stp_xml_init();
1905
1906  doc = stp_mxmlLoadFile(NULL, fp, STP_MXML_NO_CALLBACK);
1907
1908  curve = xml_doc_get_curve(doc);
1909
1910  if (doc)
1911    stp_mxmlDelete(doc);
1912
1913  stp_xml_exit();
1914  return curve;
1915
1916}
1917
1918stp_curve_t *
1919stp_curve_create_from_string(const char* string)
1920{
1921  stp_curve_t *curve = NULL;
1922  stp_mxml_node_t *doc;
1923  stp_deprintf(STP_DBG_XML,
1924	       "stp_curve_create_from_string: reading '%s'...\n", string);
1925  stp_xml_init();
1926
1927  doc = stp_mxmlLoadString(NULL, string, STP_MXML_NO_CALLBACK);
1928
1929  curve = xml_doc_get_curve(doc);
1930
1931  if (doc)
1932    stp_mxmlDelete(doc);
1933
1934  stp_xml_exit();
1935  return curve;
1936}
1937