Logo Search packages:      
Sourcecode: sagcad version File versions  Download package

spline.c

/* ====================================================================
 * ===  Copyright (C) 1998-2007 Yutaka Sagiya. All rights reserved. ===
 * ====================================================================
 * 
 *    Project              : SagCAD
 *    Description          : CAD/CAM
 *    Source               : spline.c
 * 
 *    ----------------------------------
 * 
 *    License              : GNU General Public License (GPL)
 *    Copyright            : (C) 1998-2007 by Yutaka Sagiya
 *    email                : kappa@a6s.highway.ne.jp
 *                         : yutaka@sagiya.com
 *    Begin                : 2003/01/06
 *    Last                 : 2007/11/08
 * ====================================================================
 */

#ifdef HAVE_CONFIG_H
#  include "config.h"
#endif

#include <math.h>
#include <stdlib.h>
#include "MemoryLeak.h"
#include "List_PolyLine.h"
#define _SPLINE_
#include "spline.h"



//スプライン補間


/*
平面上の任意のN点を通る閉曲線を求める。
まず x[0..N-1],y[0..N-1]にN点の座標を与え maketable2()を1度だけ呼び出す。
(条件x{0}<x{1}<...<x{N-1}は不要)。 次に 0≦t≦1の範囲で媒介変数 を動かしながら
spline2(t,&u,&v,p,x,y,a,b)を何度も呼び出し (u,v)をプロットすれば (x{0},y{0})
から出発して (x{1},y{1}),(x{2},y{2}),...を順に通り (x{N-1},y{N-1})に至る曲線が
描ける。
*/
/* 開曲線 */

//#define  N  5
//double x[N] = { 1, 2, 4, 6, 5 },
//           y[N] = { 1, 3, 4, 3, 1 },
//           p[N], a[N], b[N];
//#define  N      5
//double x[N] = { 0, 10, 20, 30, 40 },
//       y[N] = { 610.66, 1227.4, 2338.1, 4244.9, 7381.2 },
//       p[N], a[N], b[N];



void open_maketable(int n, double x[], double y[], double z[], double h[], double d[])
{
      int i;
      double t;
//    static double h[n], d[n];

      z[0] = 0;  z[n - 1] = 0;  /* 両端点での y''(x) / 6 */
      for (i = 0; i < n - 1; i++) {
            h[i    ] =  x[i + 1] - x[i];
            d[i + 1] = (y[i + 1] - y[i]) / h[i];
      }
      z[1] = d[2] - d[1] - h[0] * z[0];
      d[1] = 2 * (x[2] - x[0]);
      for (i = 1; i < n - 2; i++) {
            t = h[i] / d[i];
            z[i + 1] = d[i + 2] - d[i + 1] - z[i] * t;
            d[i + 1] = 2 * (x[i + 2] - x[i]) - h[i] * t;
      }
      z[n - 2] -= h[n - 2] * z[n - 1];
      for (i = n - 2; i > 0; i--)
            z[i] = (z[i] - h[i] * z[i + 1]) / d[i];
}



double open_spline(int n, double t, double x[], double y[], double z[])
{
      int i, j, k;
      double d, h;

      i = 0;      j = n - 1;
      while (i < j) {
            k = (i + j) / 2;
            if (x[k] < t) i = k + 1;  else j = k;
      }
      if (i > 0) i--;
      h = x[i + 1] - x[i];  d = t - x[i];
      return (((z[i + 1] - z[i]) * d / h + z[i] * 3) * d
            + ((y[i + 1] - y[i]) / h
            - (z[i] * 2 + z[i + 1]) * h)) * d + y[i];
}



void open_maketable2(int n, double p[], double x[], double y[], double a[], double b[], double h[], double d[])
{
      int i;
      double t1, t2;

      p[0] = 0;
      for (i = 1; i < n; i++) {
            t1 = x[i] - x[i - 1];
            t2 = y[i] - y[i - 1];
            p[i] = p[i - 1] + sqrt(t1 * t1 + t2 * t2);
      }
      for (i = 1 ; i < n ; i++) p[i] /= p[n - 1];
      open_maketable(n, p, x, a, h, d);
      open_maketable(n, p, y, b, h, d);
}



void open_spline2(int n, double t, double *px, double *py, double p[], double x[], double y[], double a[], double b[])
{
      *px = open_spline(n, t, p, x, a);
      *py = open_spline(n, t, p, y, b);
}



/* open spline */
#include <stdio.h>
#include <stdlib.h>
int open_spline_main(int n, VERTEX *vertex, long color, int split)
{
      int i;
      double pitch;
      static double old_x = 0, old_y = 0;

      double u, v;
      double *x, *y, *p, *a, *b;
      double *h, *d;


      x = (double *) xmalloc( n * sizeof(double) );
      y = (double *) xmalloc( n * sizeof(double) );
      p = (double *) xmalloc( n * sizeof(double) );
      a = (double *) xmalloc( n * sizeof(double) );
      b = (double *) xmalloc( n * sizeof(double) );

      h = (double *) xmalloc( n * sizeof(double) );
      d = (double *) xmalloc( n * sizeof(double) );


      for (i = 0 ; i < n ; i++) {
            x[i] = vertex[i].x;
            y[i] = vertex[i].y;
            p[i] = 0;
            a[i] = 0;
            b[i] = 0;

            h[i] = 0;
            d[i] = 0;
      }


      open_maketable2(n, p, x, y, a, b, h, d);

      pitch = 1.0/(double)split;
      for (i = 0; i <= split ; i++) {
            open_spline2(n, pitch * (double)i, &u, &v, p, x, y, a, b);
            if (i == 0) {
                  // pset
            }
            else {
                  //LineDraw(drawing_area, old_x, old_y, u, v, 1, color);
            }
            old_x = u;
            old_y = v;
//          printf("%5.2f %5.2f\n", u, v);
      }


      xfree(x);
      xfree(y);
      xfree(p);
      xfree(a);
      xfree(b);

      xfree(h);
      xfree(d);

      return EXIT_SUCCESS;
}





//スプライン補間(閉曲線)


/*
平面上のn点を通る閉曲線を求める。
このプログラムは pspline2() と maketable() と spline() および spline2.c の
maketable2() と spline2() を合わせ,配列の添字の上限を n-1 から n に変え,
maketable2() 中の <n (2か所) を ≦nにし,p[n-1] を p[n]にしたものである。
*/
/* 閉曲線 */

//#define  n      5
//double x[n + 1] = { 1, 2, 4, 6, 5, 1 },
//       y[n + 1] = { 1, 3, 4, 3, 1, 1 },
//       p[n + 1], a[n + 1], b[n + 1];

void close_maketable(int n, double x[], double y[], double z[], double h[], double d[], double w[])
{
      int i;
      double t;
//    static double h[n + 1], d[n + 1], w[n + 1];

      for (i = 0; i < n; i++) {
            h[i] = x[i + 1] - x[i];
            w[i] = (y[i + 1] - y[i]) / h[i];
      }
      w[n] = w[0];
      for (i = 1; i < n; i++) d[i] = 2 * (x[i + 1] - x[i - 1]);
      d[n] = 2 * (h[n - 1] + h[0]);
      for (i = 1; i <= n; i++) z[i] = w[i] - w[i - 1];
      w[1] = h[0];  w[n - 1] = h[n - 1];  w[n] = d[n];
      for (i = 2; i < n - 1; i++) w[i] = 0;
      for (i = 1; i < n; i++) {
            t = h[i] / d[i];
            z[i + 1] = z[i + 1] - z[i] * t;
            d[i + 1] = d[i + 1] - h[i] * t;
            w[i + 1] = w[i + 1] - w[i] * t;
      }
      w[0] = w[n];  z[0] = z[n];
      for (i = n - 2; i >= 0; i--) {
            t = h[i] / d[i + 1];
            z[i] = z[i] - z[i + 1] * t;
            w[i] = w[i] - w[i + 1] * t;
      }
      t = z[0] / w[0];  z[0] = t;  z[n] = t;
      for (i = 1; i < n; i++)
            z[i] = (z[i] - w[i] * t) / d[i];
}



double close_spline(int n, double t, double x[], double y[], double z[])
{
      int i, j, k;
      double d, h, period;

      period = x[n] - x[0];
      while (t > x[n]) t -= period;
      while (t < x[0]) t += period;
      i = 0;      j = n;
      while (i < j) {
            k = (i + j) / 2;
            if (x[k] < t) i = k + 1;  else j = k;
      }
      if (i > 0) i--;
      h = x[i + 1] - x[i];
      d = t - x[i];
      return (((z[i + 1] - z[i]) * d / h + z[i] * 3) * d
            + ((y[i + 1] - y[i]) / h
            - (z[i] * 2 + z[i + 1]) * h)) * d + y[i];
}



void close_maketable2(int n, double p[], double x[], double y[], double a[], double b[], double h[], double d[], double w[])
{
      int i;
      double t1, t2;

      p[0] = 0;
      for (i = 1; i <= n; i++) {
            t1 = x[i] - x[i - 1];
            t2 = y[i] - y[i - 1];
            p[i] = p[i - 1] + sqrt(t1 * t1 + t2 * t2);
      }
      for (i = 1; i <= n; i++) p[i] /= p[n];
      close_maketable(n, p, x, a, h, d, w);
      close_maketable(n, p, y, b, h, d, w);
}



void close_spline2(int n, double t, double *px, double *py, double p[], double x[], double      y[], double a[], double b[])
{
      *px = close_spline(n, t, p, x, a);
      *py = close_spline(n, t, p, y, b);
}





/* close spline */
//#include <stdio.h>
//#include <stdlib.h>
int close_spline_main(int n, VERTEX *vertex, long color, int split)
{
      int i;
      double pitch;
      static double old_x = 0, old_y = 0;

      double u, v;
      double *x, *y, *p, *a, *b;
      double *h, *d, *w;


      x = (double *) xmalloc( (n+1) * sizeof(double) );
      y = (double *) xmalloc( (n+1) * sizeof(double) );
      p = (double *) xmalloc( (n+1) * sizeof(double) );
      a = (double *) xmalloc( (n+1) * sizeof(double) );
      b = (double *) xmalloc( (n+1) * sizeof(double) );

      h = (double *) xmalloc( (n+1) * sizeof(double) );
      d = (double *) xmalloc( (n+1) * sizeof(double) );
      w = (double *) xmalloc( (n+1) * sizeof(double) );

      for (i = 0 ; i < (n+1) ; i++) {
            x[i] = vertex[i].x;
            y[i] = vertex[i].y;
            p[i] = 0;
            a[i] = 0;
            b[i] = 0;

            h[i] = 0;
            d[i] = 0;
            w[i] = 0;
      }


      close_maketable2(n, p, x, y, a, b, h, d, w);
      pitch = 1.0/(double)split;
      for (i = 0; i <= split ; i++) {
            close_spline2(n, pitch * i, &u, &v, p, x, y, a, b);
            if (i == 0) {
                  // pset
            }
            else {
//                LineDraw(drawing_area, old_x, old_y, u, v, 1, color);
            }
            old_x = u;
            old_y = v;
//          printf("%5.2f %5.2f\n", u, v);
      }


      xfree(x);
      xfree(y);
      xfree(p);
      xfree(a);
      xfree(b);

      xfree(h);
      xfree(d);
      xfree(w);

      return EXIT_SUCCESS;
}





void close_spline_test(void)
{
int n = 5;
int split = 60;
double x1[6] = { 1, 2, 4, 6, 5, 1 },
         y1[6] = { 1, 3, 4, 3, 1, 1 };

      int i;
      double pitch;
      static double old_x = 0, old_y = 0;

      double u, v;
      double *x, *y, *p, *a, *b;
      double *h, *d, *w;


      x = (double *) xmalloc( (n+1) * sizeof(double) );
      y = (double *) xmalloc( (n+1) * sizeof(double) );
      p = (double *) xmalloc( (n+1) * sizeof(double) );
      a = (double *) xmalloc( (n+1) * sizeof(double) );
      b = (double *) xmalloc( (n+1) * sizeof(double) );

      h = (double *) xmalloc( (n+1) * sizeof(double) );
      d = (double *) xmalloc( (n+1) * sizeof(double) );
      w = (double *) xmalloc( (n+1) * sizeof(double) );

      for (i = 0 ; i < (n+1) ; i++) {
            x[i] = x1[i];
            y[i] = y1[i];
            p[i] = 0;
            a[i] = 0;
            b[i] = 0;

            h[i] = 0;
            d[i] = 0;
            w[i] = 0;
      }


      close_maketable2(n, p, x, y, a, b, h, d, w);
      pitch = 1.0/(double)split;
      for (i = 0; i <= split ; i++) {
            close_spline2(n, pitch * i, &u, &v, p, x, y, a, b);
            if (i == 0) {
                  // pset
            }
            else {
//                LineDraw(drawing_area, old_x, old_y, u, v, 1, 0xffffff);
            }
            old_x = u;
            old_y = v;
//          printf("%5.2f %5.2f\n", u, v);
      }

      xfree(x);
      xfree(y);
      xfree(p);
      xfree(a);
      xfree(b);

      xfree(h);
      xfree(d);
      xfree(w);

      return;
}





/* ====================================================================
 * ===  Copyright (C) 1998-2007 Yutaka Sagiya. All rights reserved. ===
 * ====================================================================
 *    Project              : SagCAD
 *    Source               : spline.c
 * ====================================================================
 */

Generated by  Doxygen 1.6.0   Back to index