C++ Code zur Berechnung der Deltas beim Forward-Differencing für Bezier-Kurven (4 Punkte) 

/* --------------------------------------------------------------------------
   3d-vector 
   for operations in 3d euclidean space 
   -------------------------------------------------------------------------- */
class vector3d
{
public:
  double x;
 
double y; 
 
double z;

public:
  // construct
  vector3d (void) { }
  vector3d (
const class vector3d &v) { 
    x=v.x; y=v.y; z=v.z; }
  vector3d (
double x0, double y0, double z0) { 
    x=x0; y=y0; z=z0; }

  // set
  virtual void set (const class vector3d &v) { 
    x=v.x; y=v.y; z=v.z; }
 
virtual void set (double x0, double y0, double z0) { 
    x=x0; y=y0; z=z0; }

  // set
  class vector3d& operator = (const class vector3d &v) {
    set(v); 
    return(*this);
  }

  // dot operations
  class vector3d& operator += (const class vector3d &v) { 
    x += v.x; 
    y += v.y; 
    z += v.z; 
    return(*this); 
  }

  // ...
};

/* --------------------------------------------------------------------------
   3d-bezier curves and surfaces 
   -------------------------------------------------------------------------- */
class bezier3d {
protected:
  // create deltas for forward differencing
  virtual void forwdiff 
      (
long steps,                // number of steps p1...p4
       const class vector3d &p1// point 1 (start-point)
       const class vector3d &z12, // point 2 (intermediate 1)
       const class vector3d &z21, // point 3 (intermediate 2)
       const class vector3d &p2// point 4 (end point)
       class vector3d &delta1,    // out: delta 1 
      
class vector3d &delta2,    // out: delta 2 
      
class vector3d &delta3)    // out: delta 3
  {
   
double step1 = 1.0 / (double)(steps-1);
   
double step2 = step1*step1;
   
double step3 = step2*step1;

   
class vector3d a (-p1.x + 3.0*z12.x - 3.0*z21.x + p2.x,
                      -p1.y + 3.0*z12.y - 3.0*z21.y + p2.y,
                      -p1.z + 3.0*z12.z - 3.0*p3.z + p2.z);
   
class vector3d b (3.0*p1.x - 6.0*z12.x + 3.0*z21.x,
                      3.0*p1.y - 6.0*z12.y + 3.0*z21.y,
                      3.0*p1.z - 6.0*z12.z + 3.0*z21.z);
   
class vector3d c (-3.0*p1.x + 3.0*z12.x,
                      -3.0*p1.y + 3.0*z12.y,
                      -3.0*p1.z + 3.0*z12.z);

    delta1.set (a.x*step3 + b.x*step2 + c.x*step1,
                a.y*step3 + b.y*step2 + c.y*step1,
                a.z*step3 + b.z*step2 + c.z*step1);
    delta2.set (6.0*a.x*step3 + 2.0*b.x*step2,
                6.0*a.y*step3 + 2.0*b.y*step2,
                6.0*a.z*step3 + 2.0*b.z*step2);
    delta3.set (6.0*a.x*step3,
                6.0*a.y*step3,
                6.0*a.z*step3);
  }

 
// ...
}


C++ Code zur Berechnung der Zwischenpunkte aus den zwei Eckpunkten und zwei Normalenvektoren

/* --------------------------------------------------------------------------
   3d-bezier-curve
   2 points and 2 normal vectors
   -------------------------------------------------------------------------- */
class bezier_curve_2p2n : public bezier_curve_4p {

public:

  // create intermediat points p1...p2: p1 p12 p21 p2
  virtual void intermediates (const class vector3d &p1, // point 1
                              const class vector3d &n1, // normal vector at p1
                              const class vector3d &p2, // point 2
                              const class vector3d &n2, // normal vector at p2
                              class vector3d &p12,      // (OUT)
                              class vector3d &p21,      // (OUT)
                              double h=1.0) const       // distance factor
  {
    // lines p2+t*n1 / p1+t*n2
    class line3d l1 (p2,n1);
    class line3d l2 (p1,n2);

    // tangential planes for p1 / p2
    class plane3d e1 (n1,p1);
    class plane3d e2 (n2,p2);

    // dist |p1-p2|
    double d12 = p1.dist (p2);
 
    // intermediate points: p1 p12 p21 p2
    p12 = e1.intersection (l1);
    p21 = e2.intersection (l2);

    class vector3d v12 (p2-p1);
    
    // set point distance to max 1/3 of |p1-p2|
    double f = 1.0/(3.0*ABS(v12.cosinus(p1-p12))) * h;
    double g = 1.0/(3.0*ABS(v12.cosinus(p2-p21))) * h;
    p12 -= p1; p12 *= MIN(1.0,(d12/p12.length())*f); p12 += p1;
    p21 -= p2; p21 *= MIN(1.0,(d12/p21.length())*g); p21 += p2;
    
    // take p1 / p2 itself when p12 / p21 is not valid
    if (!p12.isvalid()) p12.set (p1); 
    if (!p21.isvalid()) p21.set (p2); 
  }
};

C++ Code zur Berechnung der Deltas beim Forward-Differencing für Bezier-Oberflächen (4x4 Punkte)

/* --------------------------------------------------------------------------
   3d-bezier-surface
   4x4 points
   -------------------------------------------------------------------------- */
class bezier_surface_4x4p {

public:

  class vector3d point  // current point

protected:

  class vector3d delta1;
  class vector3d delta2;
  class vector3d delta3;
  class vector3d delta[4][4];

protected:

  // calculate forward differences
  virtual void forwdiff (long xsteps, long ysteps
const class vector3d p[4][4]) 
  {
    double x  = 1.0 / (double)(xsteps-1); 
    double x2 = x*x, x3 = x2*x;
    
    double y = 1.0 / (double)(ysteps-1);
    double y2 = y*y, y3 = y2*y;
    
    double s[4][4] = {
      {  1.00,  0.00,  0.00,  0.00 },
      { ( -1.0*x3 + 3.0*x2 - 3.0*x), (  3.0*x3 - 6.0*x2 + 3.0*x),
        ( -3.0*x3 + 3.0*x2),         (  1.0*x3) },
      { ( -6.0*x3 + 6.0*x2),         ( 18.0*x3 - 12.0*x2), 
        (-18.0*x3 + 6.0*x2),         (  6.0*x3) },
      { ( -6.0*x3), ( 18.0*x3), (-18.0*x3), (  6.0*x3) }
    };
    double t[4][4] = {
      { 1.00, (-1.0*y3 + 3.0*y2 - 3.0*y), (-6.0*y3 +  6.0*y2), ( -6.0*y3) },
      { 0.00, ( 3.0*y3 - 6.0*y2 + 3.0*y), (18.0*y3 - 12.0*y2), ( 18.0*y3) },
      { 0.00, (-3.0*y3 + 3.0*y2),         (-18.0*y3 + 6.0*y2), (-18.0*y3) },
      { 0.00, (1.0*y3),                   ( 6.0*y3),           (  6.0*y3) }
    };

    class vector3d tdelta[4][4];
    matmult (p, t,     tdelta);
    matmult (s, tdelta, delta);
    
    point  = delta[0][0];
    delta1 = delta[0][1];
    delta2 = delta[0][2]; 
    delta3 = delta[0][3];    
  }

  // matrix multiplication c = a * b
  virtual void matmult (double a[4][4], const class vector3d b[4][4], 
class vector3d c[4][4]) {
    for (long i,y=0; y<4; y++)
      for (long x=0; x<4; x++) 
  for (i=0, c[x][y]=0; i<4; i++) { c[x][y] += (b[i][y]*a[x][i]); }
  }
  virtual void matmult (const class vector3d a[4][4], double b[4][4], 
class vector3d c[4][4]) {
    for (long i,y=0; y<4; y++)
      for (long x=0; x<4; x++) 
  for (i=0, c[x][y]=0; i<4; i++) { c[x][y] += (a[x][i]*b[i][y]); }
  }

public :

  // cunstruct / set
  bezier_surface_4x4p (void) { }
  bezier_surface_4x4p (const bezier_surface_4x4p &s) { set (s); }
  bezier_surface_4x4p (const class vector3d p[4][4],
long xsteps, long ysteps) { set (p,xsteps,ysteps); }
  virtual void set (const bezier_surface_4x4p &s) {
    point = s.point; delta1 = s.delta1; delta2 = s.delta2; delta3 = s.delta3; 
    for (long y=0; y<4; y++
      for (long x=0; x<4; x++) delta[x][y] = s.delta[x][y]; }
  virtual void set (const class vector3d p[4][4],
long xsteps, long ysteps) { forwdiff (xsteps,ysteps,p); }

  // next row / column (rows first)
  virtual void nextrow (void) {
    point += delta1; delta1 += delta2; delta2 += delta3; }
  virtual void nextcol (void) {
    for (long y=0; y<4; y++) 
      for (long x=0; x<3; x++) delta[x][y] += delta[x+1][y];
    point  = delta[0][0];
    delta1 = delta[0][1];
    delta2 = delta[0][2]; 
    delta3 = delta[0][3];
  }
};


 

C++ Code zur Berechnung der Hilfspunkte einer Bezier-Fläche

/* --------------------------------------------------------------------------
   3d-bezier-surface
   4 points + 4 normal vectors
   -------------------------------------------------------------------------- */

class bezier_surface_4p4n : public bezier_surface_4x4p {

protected:


  // border intermediates

  // n : normals for corners p[0][0], p[3][0], p[3][3], p[0][3]

  // h : intermediate point distance factor (<1:flat, >1:round)

  // p : pre: corners filled in (IN/OUT)

  virtual void outer (const class vector3d n[4], doubleh[4],

                      class vector3d p[4][4]) const

  {

    // ... border

    class bezier_curve_2p2n bc;

    bc.intermediates (p[0][0],n[0], p[3][0],n[1], p[1][0],p[2][0], h[0]);

    bc.intermediates (p[0][3],n[3], p[3][3],n[2], p[1][3],p[2][3], h[2]);

    bc.intermediates (p[0][0],n[0], p[0][3],n[3], p[0][1],p[0][2], h[3]);

    bc.intermediates (p[3][0],n[1], p[3][3],n[2], p[3][1],p[3][2], h[1]);

  }


   // inner intermediates
  // p : pre: corners / border points filled in (IN/OUT)

  virtual void inner (class vector3d p[4][4]) const

  {

    p[1][1] = inner (p[0][3],p[0][1],p[0][0],p[1][0],p[3][0],

                     p[3][1],p[3][3],p[1][3]); 

    p[2][1] = inner (p[0][0],p[2][0],p[3][0],p[3][1],p[3][3],

                     p[2][3],p[0][3],p[0][1]);

    p[2][2] = inner (p[3][0],p[3][2],p[3][3],p[2][3],p[0][3],

                     p[0][2],p[0][0],p[2][0]); 

    p[1][2] = inner (p[3][3],p[1][3],p[0][3],p[0][2],p[0][0],

                     p[1][0],p[3][0],p[3][2]);

  }


  // inner intermediate point for point 2

  //

  // p1 ... p12 p2

  //        res p23

  //            ...

  //            p3

  virtual class vector3d inner (const class vector3d &p1,

                                const class vector3d &p12,

                                const class vector3d &p2

                                const class vector3d &p23,

                                const class vector3d &p3,

                                const class vector3d &p34,

                                const class vector3d &p4,

                                const class vector3d &p14) const

  {

    // use p4 and p14 instead of p1 / p3 if p2==p1 / p2==p3

    const class vector3d *v1 = &p1, *v12 = &p12, *v3 = &p3, *v23 = &p23;

    if (p1==p2) { v1 = &p4; v12 = &p14; }

    if (p2==p3) { v3 = &p4; v23 = &p34; }


    // lines p1 - p2 / p2 -p3

    class line3d l12 (p2,p2-*v1);

    class line3d l23 (p2,p2-*v3);


    // projections border-intermediate-points on to lines l12 / l23

    class vector3d i12 (l12.projection (*v12));

    class vector3d i23 (l23.projection (*v23));


    // tangential planes at p12 / p23

    class plane3d e12 (i12-*v12,*v12);

    class plane3d e23 (i23-*v23,*v23);


    // intersection line between planes e12 / e23

    class line3d l (e12.intersection (e23));


    // extend p2 - v12 and p2 - v23

    class vector3d c12 (*v12-p2);

    class vector3d c23 (*v23-p2);

    double f = 2.0 - c12.cosinus(c23);

    class vector3d x12 (p2 + (c12*f)); 

    class vector3d x23 (p2 + (c23*f));


    // result

    class vector3d res ((l.projection (x12)+l.projection (x23))/2.0);

    if (!res.isvalid()) res = (*v12+*v23-p2);

    return (res);

  }


public:


  // create intermediate points 

  // construct / set

  // p      : corner points   p[0] p[1]

  //                          p[3] p[2]

  // n      : normal vectors at corner points

  // h      : intermediate point distance factor (<1:flat, >1:round)

  virtual void intermediates (const class vector3d p[4],

                              const class vector3d n[4],

                              class vector3d p16[4][4]) const

  {

    double h[4]={1.0,1.0,1.0,1.0}; 

    intermediates (p,n,p16,h); 

  }

  virtual void intermediates (const class vector3d p[4],

                              const class vector3d n[4],

                              class vector3d p16[4][4],

                              double h[4]) const

  {

    // corners

    p16[0][0] = p[0]; 

    p16[3][0] = p[1]; 

    p16[0][3] = p[3]; 

    p16[3][3] = p[2]; 

    // intermediates ...

    // ... border

    outer (n,h,p16);

    // ... inner

    inner (p16);

  }

  // ...

};