summaryrefslogtreecommitdiffstats
path: root/tkblt/generic/tkbltGrElemLineSpline.C
diff options
context:
space:
mode:
Diffstat (limited to 'tkblt/generic/tkbltGrElemLineSpline.C')
-rw-r--r--tkblt/generic/tkbltGrElemLineSpline.C1086
1 files changed, 0 insertions, 1086 deletions
diff --git a/tkblt/generic/tkbltGrElemLineSpline.C b/tkblt/generic/tkbltGrElemLineSpline.C
deleted file mode 100644
index 9224d53..0000000
--- a/tkblt/generic/tkbltGrElemLineSpline.C
+++ /dev/null
@@ -1,1086 +0,0 @@
-/*
- * Smithsonian Astrophysical Observatory, Cambridge, MA, USA
- * This code has been modified under the terms listed below and is made
- * available under the same terms.
- */
-
-/*
- * Copyright 2009 George A Howlett.
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject to
- * the following conditions:
- *
- * The above copyright notice and this permission notice shall be
- * included in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
- * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
- * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
- * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
- * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
- *
- */
-
-#include <float.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include <cmath>
-
-#include "tkbltGrElemLine.h"
-
-using namespace Blt;
-
-typedef double TriDiagonalMatrix[3];
-typedef struct {
- double b, c, d;
-} Cubic2D;
-
-typedef struct {
- double b, c, d, e, f;
-} Quint2D;
-
-// Quadratic spline parameters
-#define E1 param[0]
-#define E2 param[1]
-#define V1 param[2]
-#define V2 param[3]
-#define W1 param[4]
-#define W2 param[5]
-#define Z1 param[6]
-#define Z2 param[7]
-#define Y1 param[8]
-#define Y2 param[9]
-
-/*
- *---------------------------------------------------------------------------
- *
- * Search --
- *
- * Conducts a binary search for a value. This routine is called
- * only if key is between x(0) and x(len - 1).
- *
- * Results:
- * Returns the index of the largest value in xtab for which
- * x[i] < key.
- *
- *---------------------------------------------------------------------------
- */
-static int Search(Point2d points[], int nPoints, double key, int *foundPtr)
-{
- int low = 0;
- int high = nPoints - 1;
-
- while (high >= low) {
- int mid = (high + low) / 2;
- if (key > points[mid].x)
- low = mid + 1;
- else if (key < points[mid].x)
- high = mid - 1;
- else {
- *foundPtr = 1;
- return mid;
- }
- }
- *foundPtr = 0;
- return low;
-}
-
-/*
- *---------------------------------------------------------------------------
- *
- * QuadChoose --
- *
- * Determines the case needed for the computation of the parame-
- * ters of the quadratic spline.
- *
- * Results:
- * Returns a case number (1-4) which controls how the parameters
- * of the quadratic spline are evaluated.
- *
- *---------------------------------------------------------------------------
- */
-static int QuadChoose(Point2d* p, Point2d* q, double m1, double m2,
- double epsilon)
-{
- // Calculate the slope of the line joining P and Q
- double slope = (q->y - p->y) / (q->x - p->x);
-
- if (slope != 0.0) {
- double prod1 = slope * m1;
- double prod2 = slope * m2;
-
- // Find the absolute values of the slopes slope, m1, and m2
- double mref = fabs(slope);
- double mref1 = fabs(m1);
- double mref2 = fabs(m2);
-
- // If the relative deviation of m1 or m2 from slope is less than
- // epsilon, then choose case 2 or case 3.
- double relerr = epsilon * mref;
- if ((fabs(slope - m1) > relerr) && (fabs(slope - m2) > relerr) &&
- (prod1 >= 0.0) && (prod2 >= 0.0)) {
- double prod = (mref - mref1) * (mref - mref2);
- if (prod < 0.0) {
- // l1, the line through (x1,y1) with slope m1, and l2,
- // the line through (x2,y2) with slope m2, intersect
- // at a point whose abscissa is between x1 and x2.
- // The abscissa becomes a knot of the spline.
- return 1;
- }
- if (mref1 > (mref * 2.0)) {
- if (mref2 <= ((2.0 - epsilon) * mref))
- return 3;
- }
- else if (mref2 <= (mref * 2.0)) {
- // Both l1 and l2 cross the line through
- // (x1+x2)/2.0,y1 and (x1+x2)/2.0,y2, which is the
- // midline of the rectangle formed by P and Q or both
- // m1 and m2 have signs different than the sign of
- // slope, or one of m1 and m2 has opposite sign from
- // slope and l1 and l2 intersect to the left of x1 or
- // to the right of x2. The point (x1+x2)/2. is a knot
- // of the spline.
- return 2;
- }
- else if (mref1 <= ((2.0 - epsilon) * mref)) {
- // In cases 3 and 4, sign(m1)=sign(m2)=sign(slope).
- // Either l1 or l2 crosses the midline, but not both.
- // Choose case 4 if mref1 is greater than
- // (2.-epsilon)*mref; otherwise, choose case 3.
- return 3;
- }
- // If neither l1 nor l2 crosses the midline, the spline
- // requires two knots between x1 and x2.
- return 4;
- }
- else {
- // The sign of at least one of the slopes m1 or m2 does not
- // agree with the sign of *slope*.
- if ((prod1 < 0.0) && (prod2 < 0.0)) {
- return 2;
- }
- else if (prod1 < 0.0) {
- if (mref2 > ((epsilon + 1.0) * mref))
- return 1;
- else
- return 2;
- }
- else if (mref1 > ((epsilon + 1.0) * mref))
- return 1;
- else
- return 2;
- }
- }
- else if ((m1 * m2) >= 0.0)
- return 2;
- else
- return 1;
-}
-
-/*
- *---------------------------------------------------------------------------
- * Computes the knots and other parameters of the spline on the
- * interval PQ.
- * On input--
- * P and Q are the coordinates of the points of interpolation.
- * m1 is the slope at P.
- * m2 is the slope at Q.
- * ncase controls the number and location of the knots.
- * On output--
- *
- * (v1,v2),(w1,w2),(z1,z2), and (e1,e2) are the coordinates of
- * the knots and other parameters of the spline on P.
- * (e1,e2) and Q are used only if ncase=4.
- *---------------------------------------------------------------------------
- */
-static void QuadCases(Point2d* p, Point2d* q, double m1, double m2,
- double param[], int which)
-{
- if ((which == 3) || (which == 4)) {
- double c1 = p->x + (q->y - p->y) / m1;
- double d1 = q->x + (p->y - q->y) / m2;
- double h1 = c1 * 2.0 - p->x;
- double j1 = d1 * 2.0 - q->x;
- double mbar1 = (q->y - p->y) / (h1 - p->x);
- double mbar2 = (p->y - q->y) / (j1 - q->x);
-
- if (which == 4) {
- // Case 4
- Y1 = (p->x + c1) / 2.0;
- V1 = (p->x + Y1) / 2.0;
- V2 = m1 * (V1 - p->x) + p->y;
- Z1 = (d1 + q->x) / 2.0;
- W1 = (q->x + Z1) / 2.0;
- W2 = m2 * (W1 - q->x) + q->y;
- double mbar3 = (W2 - V2) / (W1 - V1);
- Y2 = mbar3 * (Y1 - V1) + V2;
- Z2 = mbar3 * (Z1 - V1) + V2;
- E1 = (Y1 + Z1) / 2.0;
- E2 = mbar3 * (E1 - V1) + V2;
- }
- else {
- // Case 3
- double k1 = (p->y - q->y + q->x * mbar2 - p->x * mbar1) / (mbar2 - mbar1);
- if (fabs(m1) > fabs(m2)) {
- Z1 = (k1 + p->x) / 2.0;
- } else {
- Z1 = (k1 + q->x) / 2.0;
- }
- V1 = (p->x + Z1) / 2.0;
- V2 = p->y + m1 * (V1 - p->x);
- W1 = (q->x + Z1) / 2.0;
- W2 = q->y + m2 * (W1 - q->x);
- Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
- }
- }
- else if (which == 2) {
- // Case 2
- Z1 = (p->x + q->x) / 2.0;
- V1 = (p->x + Z1) / 2.0;
- V2 = p->y + m1 * (V1 - p->x);
- W1 = (Z1 + q->x) / 2.0;
- W2 = q->y + m2 * (W1 - q->x);
- Z2 = (V2 + W2) / 2.0;
- }
- else {
- // Case 1
- Z1 = (p->y - q->y + m2 * q->x - m1 * p->x) / (m2 - m1);
- double ztwo = p->y + m1 * (Z1 - p->x);
- V1 = (p->x + Z1) / 2.0;
- V2 = (p->y + ztwo) / 2.0;
- W1 = (Z1 + q->x) / 2.0;
- W2 = (ztwo + q->y) / 2.0;
- Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
- }
-}
-
-static int QuadSelect(Point2d* p, Point2d* q, double m1, double m2,
- double epsilon, double param[])
-{
- int ncase = QuadChoose(p, q, m1, m2, epsilon);
- QuadCases(p, q, m1, m2, param, ncase);
- return ncase;
-}
-
-static double QuadGetImage(double p1, double p2, double p3, double x1,
- double x2, double x3)
-{
- double A = x1 - x2;
- double B = x2 - x3;
- double C = x1 - x3;
-
- double y = (p1 * (A * A) + p2 * 2.0 * B * A + p3 * (B * B)) / (C * C);
- return y;
-}
-
-/*
- *---------------------------------------------------------------------------
- * Finds the image of a point in x.
- * On input
- * x Contains the value at which the spline is evaluated.
- * leftX, leftY
- * Coordinates of the left-hand data point used in the
- * evaluation of x values.
- * rightX, rightY
- * Coordinates of the right-hand data point used in the
- * evaluation of x values.
- * Z1, Z2, Y1, Y2, E2, W2, V2
- * Parameters of the spline.
- * ncase Controls the evaluation of the spline by indicating
- * whether one or two knots were placed in the interval
- * (xtabs,xtabs1).
- *---------------------------------------------------------------------------
- */
-static void QuadSpline(Point2d* intp, Point2d* left, Point2d* right,
- double param[], int ncase)
-
-{
- double y;
-
- if (ncase == 4) {
- // Case 4: More than one knot was placed in the interval.
- // Determine the location of data point relative to the 1st knot.
- if (Y1 > intp->x)
- y = QuadGetImage(left->y, V2, Y2, Y1, intp->x, left->x);
- else if (Y1 < intp->x) {
- // Determine the location of the data point relative to the 2nd knot.
- if (Z1 > intp->x)
- y = QuadGetImage(Y2, E2, Z2, Z1, intp->x, Y1);
- else if (Z1 < intp->x)
- y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
- else
- y = Z2;
- }
- else
- y = Y2;
- }
- else {
- // Cases 1, 2, or 3:
- // Determine the location of the data point relative to the knot.
- if (Z1 < intp->x)
- y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
- else if (Z1 > intp->x)
- y = QuadGetImage(left->y, V2, Z2, Z1, intp->x, left->x);
- else
- y = Z2;
- }
-
- intp->y = y;
-}
-
-/*
- *---------------------------------------------------------------------------
- * Calculates the derivative at each of the data points. The
- * slopes computed will insure that an osculatory quadratic
- * spline will have one additional knot between two adjacent
- * points of interpolation. Convexity and monotonicity are
- * preserved wherever these conditions are compatible with the
- * data.
- *---------------------------------------------------------------------------
- */
-static void QuadSlopes(Point2d *points, double *m, int nPoints)
-{
- double m1s =0;
- double m2s =0;
- double m1 =0;
- double m2 =0;
- int i, n, l;
- for (l = 0, i = 1, n = 2; i < (nPoints - 1); l++, i++, n++) {
- // Calculate the slopes of the two lines joining three
- // consecutive data points.
- double ydif1 = points[i].y - points[l].y;
- double ydif2 = points[n].y - points[i].y;
- m1 = ydif1 / (points[i].x - points[l].x);
- m2 = ydif2 / (points[n].x - points[i].x);
- if (i == 1) {
- // Save slopes of starting point
- m1s = m1;
- m2s = m2;
- }
- // If one of the preceding slopes is zero or if they have opposite
- // sign, assign the value zero to the derivative at the middle point.
- if ((m1 == 0.0) || (m2 == 0.0) || ((m1 * m2) <= 0.0))
- m[i] = 0.0;
- else if (fabs(m1) > fabs(m2)) {
- // Calculate the slope by extending the line with slope m1.
- double xbar = ydif2 / m1 + points[i].x;
- double xhat = (xbar + points[n].x) / 2.0;
- m[i] = ydif2 / (xhat - points[i].x);
- }
- else {
- // Calculate the slope by extending the line with slope m2.
- double xbar = -ydif1 / m2 + points[i].x;
- double xhat = (points[l].x + xbar) / 2.0;
- m[i] = ydif1 / (points[i].x - xhat);
- }
- }
-
- // Calculate the slope at the last point, x(n).
- i = nPoints - 2;
- n = nPoints - 1;
- if ((m1 * m2) < 0.0)
- m[n] = m2 * 2.0;
- else {
- double xmid = (points[i].x + points[n].x) / 2.0;
- double yxmid = m[i] * (xmid - points[i].x) + points[i].y;
- m[n] = (points[n].y - yxmid) / (points[n].x - xmid);
- if ((m[n] * m2) < 0.0)
- m[n] = 0.0;
- }
-
- // Calculate the slope at the first point, x(0).
- if ((m1s * m2s) < 0.0)
- m[0] = m1s * 2.0;
- else {
- double xmid = (points[0].x + points[1].x) / 2.0;
- double yxmid = m[1] * (xmid - points[1].x) + points[1].y;
- m[0] = (yxmid - points[0].y) / (xmid - points[0].x);
- if ((m[0] * m1s) < 0.0)
- m[0] = 0.0;
- }
-}
-
-/*
- *---------------------------------------------------------------------------
- *
- * QuadEval --
- *
- * QuadEval controls the evaluation of an osculatory quadratic
- * spline. The user may provide his own slopes at the points of
- * interpolation or use the subroutine 'QuadSlopes' to calculate
- * slopes which are consistent with the shape of the data.
- *
- * ON INPUT--
- * intpPts must be a nondecreasing vector of points at which the
- * spline will be evaluated.
- * origPts contains the abscissas of the data points to be
- * interpolated. xtab must be increasing.
- * y contains the ordinates of the data points to be
- * interpolated.
- * m contains the slope of the spline at each point of
- * interpolation.
- * nPoints number of data points (dimension of xtab and y).
- * numEval is the number of points of evaluation (dimension of
- * xval and yval).
- * epsilon is a relative error tolerance used in subroutine
- * 'QuadChoose' to distinguish the situation m(i) or
- * m(i+1) is relatively close to the slope or twice
- * the slope of the linear segment between xtab(i) and
- * xtab(i+1). If this situation occurs, roundoff may
- * cause a change in convexity or monotonicity of the
- * resulting spline and a change in the case number
- * provided by 'QuadChoose'. If epsilon is not equal to zero,
- * then epsilon should be greater than or equal to machine
- * epsilon.
- * ON OUTPUT--
- * yval contains the images of the points in xval.
- * err is one of the following error codes:
- * 0 - QuadEval ran normally.
- * 1 - xval(i) is less than xtab(1) for at least one
- * i or xval(i) is greater than xtab(num) for at
- * least one i. QuadEval will extrapolate to provide
- * function values for these abscissas.
- * 2 - xval(i+1) < xval(i) for some i.
- *
- *
- * QuadEval calls the following subroutines or functions:
- * Search
- * QuadCases
- * QuadChoose
- * QuadSpline
- *---------------------------------------------------------------------------
- */
-static int QuadEval(Point2d origPts[], int nOrigPts, Point2d intpPts[],
- int nIntpPts, double *m, double epsilon)
-{
- double param[10];
-
- // Initialize indices and set error result
- int error = 0;
- int l = nOrigPts - 1;
- int p = l - 1;
- int ncase = 1;
-
- // Determine if abscissas of new vector are non-decreasing.
- for (int jj=1; jj<nIntpPts; jj++) {
- if (intpPts[jj].x < intpPts[jj - 1].x)
- return 2;
- }
- // Determine if any of the points in xval are LESS than the
- // abscissa of the first data point.
- int start;
- for (start = 0; start < nIntpPts; start++) {
- if (intpPts[start].x >= origPts[0].x)
- break;
- }
- // Determine if any of the points in xval are GREATER than the
- // abscissa of the l data point.
- int end;
- for (end = nIntpPts - 1; end >= 0; end--) {
- if (intpPts[end].x <= origPts[l].x)
- break;
- }
-
- if (start > 0) {
- // Set error value to indicate that extrapolation has occurred
- error = 1;
-
- // Calculate the images of points of evaluation whose abscissas
- // are less than the abscissa of the first data point.
- ncase = QuadSelect(origPts, origPts + 1, m[0], m[1], epsilon, param);
- for (int jj=0; jj<(start - 1); jj++)
- QuadSpline(intpPts + jj, origPts, origPts + 1, param, ncase);
- if (nIntpPts == 1)
- return error;
- }
- int ii;
- int nn;
- if ((nIntpPts == 1) && (end != (nIntpPts - 1)))
- goto noExtrapolation;
-
- // Search locates the interval in which the first in-range
- // point of evaluation lies.
- int found;
- ii = Search(origPts, nOrigPts, intpPts[start].x, &found);
-
- nn = ii + 1;
- if (nn >= nOrigPts) {
- nn = nOrigPts - 1;
- ii = nOrigPts - 2;
- }
- /*
- * If the first in-range point of evaluation is equal to one
- * of the data points, assign the appropriate value from y.
- * Continue until a point of evaluation is found which is not
- * equal to a data point.
- */
- if (found) {
- do {
- intpPts[start].y = origPts[ii].y;
- start++;
- if (start >= nIntpPts) {
- return error;
- }
- } while (intpPts[start - 1].x == intpPts[start].x);
-
- for (;;) {
- if (intpPts[start].x < origPts[nn].x) {
- break; /* Break out of for-loop */
- }
- if (intpPts[start].x == origPts[nn].x) {
- do {
- intpPts[start].y = origPts[nn].y;
- start++;
- if (start >= nIntpPts) {
- return error;
- }
- } while (intpPts[start].x == intpPts[start - 1].x);
- }
- ii++;
- nn++;
- }
- }
- /*
- * Calculate the images of all the points which lie within
- * range of the data.
- */
- if ((ii > 0) || (error != 1))
- ncase = QuadSelect(origPts+ii, origPts+nn, m[ii], m[nn], epsilon, param);
-
- for (int jj=start; jj<=end; jj++) {
- // If xx(j) - x(n) is negative, do not recalculate
- // the parameters for this section of the spline since
- // they are already known.
- if (intpPts[jj].x == origPts[nn].x) {
- intpPts[jj].y = origPts[nn].y;
- continue;
- }
- else if (intpPts[jj].x > origPts[nn].x) {
- double delta;
-
- // Determine that the routine is in the correct part of the spline
- do {
- ii++;
- nn++;
- delta = intpPts[jj].x - origPts[nn].x;
- } while (delta > 0.0);
-
- if (delta < 0.0)
- ncase = QuadSelect(origPts+ii, origPts+nn, m[ii], m[nn],
- epsilon, param);
- else if (delta == 0.0) {
- intpPts[jj].y = origPts[nn].y;
- continue;
- }
- }
- QuadSpline(intpPts+jj, origPts+ii, origPts+nn, param, ncase);
- }
-
- if (end == (nIntpPts - 1))
- return error;
-
- if ((nn == l) && (intpPts[end].x != origPts[l].x))
- goto noExtrapolation;
-
- // Set error value to indicate that extrapolation has occurred
- error = 1;
- ncase = QuadSelect(origPts + p, origPts + l, m[p], m[l], epsilon, param);
-
- noExtrapolation:
- // Calculate the images of the points of evaluation whose
- // abscissas are greater than the abscissa of the last data point.
- for (int jj=(end + 1); jj<nIntpPts; jj++)
- QuadSpline(intpPts + jj, origPts + p, origPts + l, param, ncase);
-
- return error;
-}
-
-/*
- *---------------------------------------------------------------------------
- *
- * Shape preserving quadratic splines
- * by D.F.Mcallister & J.A.Roulier
- * Coded by S.L.Dodd & M.Roulier
- * N.C.State University
- *
- *---------------------------------------------------------------------------
- */
-/*
- * Driver routine for quadratic spline package
- * On input--
- * X,Y Contain n-long arrays of data (x is increasing)
- * XM Contains m-long array of x values (increasing)
- * eps Relative error tolerance
- * n Number of input data points
- * m Number of output data points
- * On output--
- * work Contains the value of the first derivative at each data point
- * ym Contains the interpolated spline value at each data point
- */
-int LineElement::quadraticSpline(Point2d *origPts, int nOrigPts,
- Point2d *intpPts, int nIntpPts)
-{
- double* work = new double[nOrigPts];
- double epsilon = 0.0;
- /* allocate space for vectors used in calculation */
- QuadSlopes(origPts, work, nOrigPts);
- int result = QuadEval(origPts, nOrigPts, intpPts, nIntpPts, work, epsilon);
- delete [] work;
- if (result > 1) {
- return 0;
- }
- return 1;
-}
-
-/*
- *---------------------------------------------------------------------------
- * Reference:
- * Numerical Analysis, R. Burden, J. Faires and A. Reynolds.
- * Prindle, Weber & Schmidt 1981 pp 112
- *---------------------------------------------------------------------------
- */
-int LineElement::naturalSpline(Point2d *origPts, int nOrigPts,
- Point2d *intpPts, int nIntpPts)
-{
- Point2d *ip, *iend;
- double x, dy, alpha;
- int isKnot;
- int i, j, n;
-
- double* dx = new double[nOrigPts];
- /* Calculate vector of differences */
- for (i = 0, j = 1; j < nOrigPts; i++, j++) {
- dx[i] = origPts[j].x - origPts[i].x;
- if (dx[i] < 0.0) {
- return 0;
- }
- }
- n = nOrigPts - 1; /* Number of intervals. */
- TriDiagonalMatrix* A = new TriDiagonalMatrix[nOrigPts];
- if (!A) {
- delete [] dx;
- return 0;
- }
- /* Vectors to solve the tridiagonal matrix */
- A[0][0] = A[n][0] = 1.0;
- A[0][1] = A[n][1] = 0.0;
- A[0][2] = A[n][2] = 0.0;
-
- /* Calculate the intermediate results */
- for (i = 0, j = 1; j < n; j++, i++) {
- alpha = 3.0 * ((origPts[j + 1].y / dx[j]) - (origPts[j].y / dx[i]) -
- (origPts[j].y / dx[j]) + (origPts[i].y / dx[i]));
- A[j][0] = 2 * (dx[j] + dx[i]) - dx[i] * A[i][1];
- A[j][1] = dx[j] / A[j][0];
- A[j][2] = (alpha - dx[i] * A[i][2]) / A[j][0];
- }
-
- Cubic2D* eq = new Cubic2D[nOrigPts];
- if (!eq) {
- delete [] A;
- delete [] dx;
- return 0;
- }
- eq[0].c = eq[n].c = 0.0;
- for (j = n, i = n - 1; i >= 0; i--, j--) {
- eq[i].c = A[i][2] - A[i][1] * eq[j].c;
- dy = origPts[i+1].y - origPts[i].y;
- eq[i].b = (dy) / dx[i] - dx[i] * (eq[j].c + 2.0 * eq[i].c) / 3.0;
- eq[i].d = (eq[j].c - eq[i].c) / (3.0 * dx[i]);
- }
- delete [] A;
- delete [] dx;
-
- /* Now calculate the new values */
- for (ip = intpPts, iend = ip + nIntpPts; ip < iend; ip++) {
- ip->y = 0.0;
- x = ip->x;
-
- /* Is it outside the interval? */
- if ((x < origPts[0].x) || (x > origPts[n].x)) {
- continue;
- }
- /* Search for the interval containing x in the point array */
- i = Search(origPts, nOrigPts, x, &isKnot);
- if (isKnot) {
- ip->y = origPts[i].y;
- } else {
- i--;
- x -= origPts[i].x;
- ip->y = origPts[i].y + x * (eq[i].b + x * (eq[i].c + x * eq[i].d));
- }
- }
- delete [] eq;
- return 1;
-}
-
-typedef struct {
- double t; /* Arc length of interval. */
- double x; /* 2nd derivative of X with respect to T */
- double y; /* 2nd derivative of Y with respect to T */
-} CubicSpline;
-
-/*
- * The following two procedures solve the special linear system which arise
- * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
- * equations can be written as (i=0,1,...,n-1):
- * m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
- * In matrix notation one gets A * x = b, where the matrix A is tridiagonal
- * with additional elements in the upper right and lower left position:
- * A[i][0] = A_{i,i-1} for i=1,2,...,n-1 and m[0][0] = A_{0,n-1} ,
- * A[i][1] = A_{i, i } for i=0,1,...,n-1
- * A[i][2] = A_{i,i+1} for i=0,1,...,n-2 and m[n-1][2] = A_{n-1,0}.
- * A should be symmetric (A[i+1][0] == A[i][2]) and positive definite.
- * The size of the system is given in n (n>=1).
- *
- * In the first procedure the Cholesky decomposition A = C^T * D * C
- * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
- * Return TRUE if decomposition exist.
- */
-static int SolveCubic1(TriDiagonalMatrix A[], int n)
-{
- int i;
- double m_ij, m_n, m_nn, d;
-
- if (n < 1) {
- return 0; /* Dimension should be at least 1 */
- }
- d = A[0][1]; /* D_{0,0} = A_{0,0} */
- if (d <= 0.0) {
- return 0; /* A (or D) should be positive definite */
- }
- m_n = A[0][0]; /* A_{0,n-1} */
- m_nn = A[n - 1][1]; /* A_{n-1,n-1} */
- for (i = 0; i < n - 2; i++) {
- m_ij = A[i][2]; /* A_{i,1} */
- A[i][2] = m_ij / d; /* C_{i,i+1} */
- A[i][0] = m_n / d; /* C_{i,n-1} */
- m_nn -= A[i][0] * m_n; /* to get C_{n-1,n-1} */
- m_n = -A[i][2] * m_n; /* to get C_{i+1,n-1} */
- d = A[i + 1][1] - A[i][2] * m_ij; /* D_{i+1,i+1} */
- if (d <= 0.0) {
- return 0; /* Elements of D should be positive */
- }
- A[i + 1][1] = d;
- }
- if (n >= 2) { /* Complete last column */
- m_n += A[n - 2][2]; /* add A_{n-2,n-1} */
- A[n - 2][0] = m_n / d; /* C_{n-2,n-1} */
- A[n - 1][1] = d = m_nn - A[n - 2][0] * m_n; /* D_{n-1,n-1} */
- if (d <= 0.0) {
- return 0;
- }
- }
- return 1;
-}
-
-/*
- * The second procedure solves the linear system, with the Cholesky
- * decomposition calculated above (in m[][]) and the right side b given
- * in x[]. The solution x overwrites the right side in x[].
- */
-static void SolveCubic2(TriDiagonalMatrix A[], CubicSpline spline[],
- int nIntervals)
-{
- int n = nIntervals - 2;
- int m = nIntervals - 1;
-
- // Division by transpose of C : b = C^{-T} * b
- double x = spline[m].x;
- double y = spline[m].y;
- for (int ii=0; ii<n; ii++) {
- spline[ii + 1].x -= A[ii][2] * spline[ii].x; /* C_{i,i+1} * x(i) */
- spline[ii + 1].y -= A[ii][2] * spline[ii].y; /* C_{i,i+1} * x(i) */
- x -= A[ii][0] * spline[ii].x; /* C_{i,n-1} * x(i) */
- y -= A[ii][0] * spline[ii].y; /* C_{i,n-1} * x(i) */
- }
- if (n >= 0) {
- // C_{n-2,n-1} * x_{n-1}
- spline[m].x = x - A[n][0] * spline[n].x;
- spline[m].y = y - A[n][0] * spline[n].y;
- }
- // Division by D: b = D^{-1} * b
- for (int ii=0; ii<nIntervals; ii++) {
- spline[ii].x /= A[ii][1];
- spline[ii].y /= A[ii][1];
- }
-
- // Division by C: b = C^{-1} * b
- x = spline[m].x;
- y = spline[m].y;
- if (n >= 0) {
- // C_{n-2,n-1} * x_{n-1}
- spline[n].x -= A[n][0] * x;
- spline[n].y -= A[n][0] * y;
- }
- for (int ii=(n - 1); ii>=0; ii--) {
- // C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1}
- spline[ii].x -= A[ii][2] * spline[ii + 1].x + A[ii][0] * x;
- spline[ii].y -= A[ii][2] * spline[ii + 1].y + A[ii][0] * y;
- }
-}
-
-/*
- * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
- * through list of points (x_i,y_i). The parameter t is calculated as the
- * length of the linear stroke. The number of points must be at least 3.
- * Note: For CLOSED_CONTOURs the first and last point must be equal.
- */
-static CubicSpline* CubicSlopes(Point2d points[], int nPoints,
- int isClosed, double unitX, double unitY)
-{
- CubicSpline *s1, *s2;
- int n, i;
- double norm, dx, dy;
-
- CubicSpline* spline = new CubicSpline[nPoints];
- if (!spline)
- return NULL;
-
- TriDiagonalMatrix *A = new TriDiagonalMatrix[nPoints];
- if (!A) {
- delete [] spline;
- return NULL;
- }
- /*
- * Calculate first differences in (dxdt2[i], y[i]) and interval lengths
- * in dist[i]:
- */
- s1 = spline;
- for (i = 0; i < nPoints - 1; i++) {
- s1->x = points[i+1].x - points[i].x;
- s1->y = points[i+1].y - points[i].y;
-
- /*
- * The Norm of a linear stroke is calculated in "normal coordinates"
- * and used as interval length:
- */
- dx = s1->x / unitX;
- dy = s1->y / unitY;
- s1->t = sqrt(dx * dx + dy * dy);
-
- s1->x /= s1->t; /* first difference, with unit norm: */
- s1->y /= s1->t; /* || (dxdt2[i], y[i]) || = 1 */
- s1++;
- }
-
- /*
- * Setup linear System: Ax = b
- */
- n = nPoints - 2; /* Without first and last point */
- if (isClosed) {
- /* First and last points must be equal for CLOSED_CONTOURs */
- spline[nPoints - 1].t = spline[0].t;
- spline[nPoints - 1].x = spline[0].x;
- spline[nPoints - 1].y = spline[0].y;
- n++; /* Add last point (= first point) */
- }
- s1 = spline, s2 = s1 + 1;
- for (i = 0; i < n; i++) {
- /* Matrix A, mainly tridiagonal with cyclic second index
- ("j = j+n mod n")
- */
- A[i][0] = s1->t; /* Off-diagonal element A_{i,i-1} */
- A[i][1] = 2.0 * (s1->t + s2->t); /* A_{i,i} */
- A[i][2] = s2->t; /* Off-diagonal element A_{i,i+1} */
-
- /* Right side b_x and b_y */
- s1->x = (s2->x - s1->x) * 6.0;
- s1->y = (s2->y - s1->y) * 6.0;
-
- /*
- * If the linear stroke shows a cusp of more than 90 degree,
- * the right side is reduced to avoid oscillations in the
- * spline:
- */
- /*
- * The Norm of a linear stroke is calculated in "normal coordinates"
- * and used as interval length:
- */
- dx = s1->x / unitX;
- dy = s1->y / unitY;
- norm = sqrt(dx * dx + dy * dy) / 8.5;
- if (norm > 1.0) {
- /* The first derivative will not be continuous */
- s1->x /= norm;
- s1->y /= norm;
- }
- s1++, s2++;
- }
-
- if (!isClosed) {
- /* Third derivative is set to zero at both ends */
- A[0][1] += A[0][0]; /* A_{0,0} */
- A[0][0] = 0.0; /* A_{0,n-1} */
- A[n-1][1] += A[n-1][2]; /* A_{n-1,n-1} */
- A[n-1][2] = 0.0; /* A_{n-1,0} */
- }
- /* Solve linear systems for dxdt2[] and y[] */
-
- if (SolveCubic1(A, n)) { /* Cholesky decomposition */
- SolveCubic2(A, spline, n); /* A * dxdt2 = b_x */
- }
- else { /* Should not happen, but who knows ... */
- delete [] A;
- delete [] spline;
- return NULL;
- }
- /* Shift all second derivatives one place right and update the ends. */
- s2 = spline + n, s1 = s2 - 1;
- for (/* empty */; s2 > spline; s2--, s1--) {
- s2->x = s1->x;
- s2->y = s1->y;
- }
- if (isClosed) {
- spline[0].x = spline[n].x;
- spline[0].y = spline[n].y;
- } else {
- /* Third derivative is 0.0 for the first and last interval. */
- spline[0].x = spline[1].x;
- spline[0].y = spline[1].y;
- spline[n + 1].x = spline[n].x;
- spline[n + 1].y = spline[n].y;
- }
- delete [] A;
- return spline;
-}
-
-// Calculate interpolated values of the spline function (defined via p_cntr
-// and the second derivatives dxdt2[] and dydt2[]). The number of tabulated
-// values is n. On an equidistant grid n_intpol values are calculated.
-static int CubicEval(Point2d *origPts, int nOrigPts, Point2d *intpPts,
- int nIntpPts, CubicSpline *spline)
-{
- double t, tSkip;
- Point2d q;
- int count;
-
- /* Sum the lengths of all the segments (intervals). */
- double tMax = 0.0;
- for (int ii=0; ii<nOrigPts - 1; ii++)
- tMax += spline[ii].t;
-
- /* Need a better way of doing this... */
-
- /* The distance between interpolated points */
- tSkip = (1. - 1e-7) * tMax / (nIntpPts - 1);
-
- t = 0.0; /* Spline parameter value. */
- q = origPts[0];
- count = 0;
-
- intpPts[count++] = q; /* First point. */
- t += tSkip;
-
- for (int ii=0, jj=1; jj<nOrigPts; ii++, jj++) {
- // Interval length
- double d = spline[ii].t;
- Point2d p = q;
- q = origPts[ii+1];
- double hx = (q.x - p.x) / d;
- double hy = (q.y - p.y) / d;
- double dx0 = (spline[jj].x + 2 * spline[ii].x) / 6.0;
- double dy0 = (spline[jj].y + 2 * spline[ii].y) / 6.0;
- double dx01 = (spline[jj].x - spline[ii].x) / (6.0 * d);
- double dy01 = (spline[jj].y - spline[ii].y) / (6.0 * d);
- while (t <= spline[ii].t) { /* t in current interval ? */
- p.x += t * (hx + (t - d) * (dx0 + t * dx01));
- p.y += t * (hy + (t - d) * (dy0 + t * dy01));
- intpPts[count++] = p;
- t += tSkip;
- }
-
- // Parameter t relative to start of next interval
- t -= spline[ii].t;
- }
-
- return count;
-}
-
-int LineElement::naturalParametricSpline(Point2d *origPts, int nOrigPts,
- Region2d *extsPtr, int isClosed,
- Point2d *intpPts, int nIntpPts)
-{
- // Generate a cubic spline curve through the points (x_i,y_i) which are
- // stored in the linked list p_cntr.
- // The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
- // parameter t is the length of the linear stroke.
-
- if (nOrigPts < 3)
- return 0;
-
- if (isClosed) {
- origPts[nOrigPts].x = origPts[0].x;
- origPts[nOrigPts].y = origPts[0].y;
- nOrigPts++;
- }
-
- // Width and height of the grid is used at unit length (2d-norm)
- double unitX = extsPtr->right - extsPtr->left;
- double unitY = extsPtr->bottom - extsPtr->top;
- if (unitX < FLT_EPSILON)
- unitX = FLT_EPSILON;
- if (unitY < FLT_EPSILON)
- unitY = FLT_EPSILON;
-
- /* Calculate parameters for cubic spline:
- * t = arc length of interval.
- * dxdt2 = second derivatives of x with respect to t,
- * dydt2 = second derivatives of y with respect to t,
- */
- CubicSpline* spline = CubicSlopes(origPts, nOrigPts, isClosed, unitX, unitY);
- if (spline == NULL)
- return 0;
-
- int result= CubicEval(origPts, nOrigPts, intpPts, nIntpPts, spline);
-
- delete [] spline;
- return result;
-}
-
-static void CatromCoeffs(Point2d* p, Point2d* a, Point2d* b,
- Point2d* c, Point2d* d)
-{
- a->x = -p[0].x + 3.0 * p[1].x - 3.0 * p[2].x + p[3].x;
- b->x = 2.0 * p[0].x - 5.0 * p[1].x + 4.0 * p[2].x - p[3].x;
- c->x = -p[0].x + p[2].x;
- d->x = 2.0 * p[1].x;
- a->y = -p[0].y + 3.0 * p[1].y - 3.0 * p[2].y + p[3].y;
- b->y = 2.0 * p[0].y - 5.0 * p[1].y + 4.0 * p[2].y - p[3].y;
- c->y = -p[0].y + p[2].y;
- d->y = 2.0 * p[1].y;
-}
-
-int LineElement::catromParametricSpline(Point2d* points, int nPoints,
- Point2d* intpPts, int nIntpPts)
-{
- // The spline is computed in screen coordinates instead of data points so
- // that we can select the abscissas of the interpolated points from each
- // pixel horizontally across the plotting area.
-
- Point2d* origPts = new Point2d[nPoints + 4];
- memcpy(origPts + 1, points, sizeof(Point2d) * nPoints);
-
- origPts[0] = origPts[1];
- origPts[nPoints + 2] = origPts[nPoints + 1] = origPts[nPoints];
-
- for (int ii=0; ii<nIntpPts; ii++) {
- int interval = (int)intpPts[ii].x;
- double t = intpPts[ii].y;
- Point2d a, b, c, d;
- CatromCoeffs(origPts + interval, &a, &b, &c, &d);
- intpPts[ii].x = (d.x + t * (c.x + t * (b.x + t * a.x))) / 2.0;
- intpPts[ii].y = (d.y + t * (c.y + t * (b.y + t * a.y))) / 2.0;
- }
-
- delete [] origPts;
- return 1;
-}
generated by cgit v0.12 at 2024-11-19 16:04:41 (GMT)