src/AppParCurves/AppParCurves.cxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #define No_Standard_RangeError
  16 #define No_Standard_OutOfRange
  17
  18
  19 #include <AppParCurves.hxx>
  20 #include <BSplCLib.hxx>
  21 #include <gp_Pnt2d.hxx>
  22 #include <gp_Vec2d.hxx>
  23 #include <math_Matrix.hxx>
  24 #include <TColStd_Array1OfReal.hxx>
  25
  26 void AppParCurves::BernsteinMatrix(const Standard_Integer NbPoles,
  27                                    const math_Vector& U,
  28                                    math_Matrix& A) {
  29
  30   Standard_Integer i, j, id;
  31   Standard_Real u0, u1, y0, y1, xs;
  32   Standard_Integer first = U.Lower(), last = U.Upper();
  33   math_Vector B(1, NbPoles-1);
  34
  35
  36   for (i = first; i <= last; i++) {
  37     B(1) = 1;
  38     u0 = U(i);
  39     u1 = 1.-u0;
  40
  41     for (id = 2; id <= NbPoles-1; id++) {
  42       y0 = B(1);
  43       y1 = u0*y0;
  44       B(1) = y0-y1;
  45       for (j = 2; j <= id-1; j++) {
  46         xs = y1;
  47         y0 = B(j);
  48         y1 = u0*y0;
  49         B(j) = y0-y1+xs;
  50       }
  51       B(id) = y1;
  52     }
  53     A(i, 1) = u1*B(1);
  54     A(i, NbPoles) = u0*B(NbPoles-1);
  55     for (j = 2; j <= NbPoles-1; j++) {
  56       A(i, j) = u1*B(j)+u0*B(j-1);
  57     }
  58   }
  59 }
  60
  61
  62 void AppParCurves::Bernstein(const Standard_Integer NbPoles,
  63                              const math_Vector& U,
  64                              math_Matrix& A,
  65                              math_Matrix& DA) {
  66
  67   Standard_Integer i, j, id, Ndeg = NbPoles-1;
  68   Standard_Real u0, u1, y0, y1, xs, bj, bj1;
  69   Standard_Integer first = U.Lower(), last = U.Upper();
  70   math_Vector B(1, NbPoles-1);
  71
  72
  73   for (i = first; i <= last; i++) {
  74     B(1) = 1;
  75     u0 = U(i);
  76     u1 = 1.-u0;
  77
  78     for (id = 2; id <= NbPoles-1; id++) {
  79       y0 = B(1);
  80       y1 = u0*y0;
  81       B(1) = y0-y1;
  82       for (j = 2; j <= id-1; j++) {
  83         xs = y1;
  84         y0 = B(j);
  85         y1 = u0*y0;
  86         B(j) = y0-y1+xs;
  87       }
  88       B(id) = y1;
  89     }
  90     DA(i, 1) = -Ndeg*B(1);
  91     DA(i, NbPoles) = Ndeg*B(NbPoles-1);
  92     A(i, 1) = u1*B(1);
  93     A(i, NbPoles) = u0*B(NbPoles-1);
  94     for (j = 2; j <= NbPoles-1; j++) {
  95       bj = B(j);  bj1 = B(j-1);
  96       DA(i,j) = Ndeg*(bj1-bj);
  97       A(i, j) = u1*bj+u0*bj1;
  98     }
  99   }
 100 }
 101
 102
 103 void AppParCurves::SecondDerivativeBernstein(const Standard_Real U,
 104                                              math_Vector&        DDA) {
 105 //  Standard_Real U1 = 1-U, Y0, Y1, Xs;
 106   Standard_Real Y0, Y1, Xs;
 107   Standard_Integer NbPoles = DDA.Length();
 108   Standard_Integer id, j, N4, deg = NbPoles-1;
 109   N4 = deg*(deg-1);
 110   math_Vector B(1, deg-1);
 111   B(1) = 1.;
 112
 113   // Cas particulier si degre = 1:
 114   if (deg == 1) {
 115     DDA(1) = 0.0;
 116     DDA(2) = 0.0;
 117   }
 118   else if (deg == 2) {
 119     DDA(1) = 2.0;
 120     DDA(2) = -4.0;
 121     DDA(3) = 2.0;
 122   }
 123   else {
 124
 125     for (id = 2; id <= deg-1; id++) {
 126       Y0 = B(1);
 127       Y1 = U*Y0;
 128       B(1) = Y0-Y1;
 129       for (j = 2; j <= id-1; j++) {
 130         Xs = Y1;
 131         Y0 = B(j);
 132         Y1 = U*Y0;
 133         B(j) = Y0 - Y1 +Xs;
 134       }
 135       B(id) = Y1;
 136     }
 137
 138     DDA(1) = N4*B(1);
 139     DDA(2) = N4*(-2*B(1) + B(2));
 140     DDA(deg) = N4*(B(deg-2) - 2*B(deg-1));
 141     DDA(deg+1) = N4*B(deg-1);
 142
 143     for(j = 2; j <= deg-2; j++) {
 144       DDA(j+1) = N4*(B(j-1)-2*B(j)+B(j+1));
 145     }
 146   }
 147 }
 148
 149
 150
 151 void AppParCurves::SplineFunction(const Standard_Integer nbpoles,
 152                                   const Standard_Integer deg,
 153                                   const math_Vector&     Parameters,
 154                                   const math_Vector&     flatknots,
 155                                   math_Matrix&           A,
 156                                   math_Matrix&           DA,
 157                                   math_IntegerVector&    index)
 158 {
 159 //  Standard_Real U, NewU, co, diff, t1, t2;
 160   Standard_Real U, NewU;
 161 //  gp_Pnt2d Pt, P0;
 162 //  gp_Vec2d V1;
 163 //  Standard_Integer i, j, k, iter, in, ik, deg1 = deg+1;
 164   Standard_Integer i, j, deg1 = deg+1;
 165 //  Standard_Integer oldkindex, kindex, theindex, ttindex;
 166   Standard_Integer oldkindex, kindex, theindex;
 167   math_Vector locpoles(1 , deg1);
 168   math_Vector locdpoles(1 , deg1);
 169   Standard_Integer firstp = Parameters.Lower(), lastp = Parameters.Upper();
 170
 171   TColStd_Array1OfReal Aflatknots(flatknots.Lower(), flatknots.Upper());
 172   for (i = flatknots.Lower(); i <= flatknots.Upper(); i++) {
 173     Aflatknots(i) = flatknots(i);
 174   }
 175
 176
 177   oldkindex = 1;
 178
 179   Standard_Integer pp, qq;
 180   Standard_Real Saved, Inverse, LocalInverse, locqq, locdqq, val;
 181
 182   for (i = firstp; i <= lastp; i++) {
 183     U = Parameters(i);
 184     NewU = U;
 185     kindex = oldkindex;
 186     BSplCLib::LocateParameter(deg, Aflatknots, U, Standard_False,
 187                               deg1, nbpoles+1, kindex, NewU);
 188
 189     oldkindex = kindex;
 190
 191     // On stocke les index:
 192     index(i) = kindex - deg-1;
 193
 194     locpoles(1) = 1.0;
 195
 196     for (qq = 2; qq <= deg; qq++) {
 197       locpoles(qq) = 0.0;
 198       for (pp = 1; pp <= qq-1; pp++) {
 199         Inverse = 1.0 / (flatknots(kindex + pp)  - flatknots(kindex - qq + pp + 1)) ;
 200         Saved = (U - flatknots(kindex - qq + pp + 1)) * Inverse * locpoles(pp);
 201         locpoles(pp) *= (flatknots(kindex + pp) - U) * Inverse ;
 202         locpoles(pp) += locpoles(qq) ;
 203         locpoles(qq) = Saved ;
 204       }
 205     }
 206
 207     qq = deg+1;
 208     for (pp = 1; pp <= deg; pp++) {
 209       locdpoles(pp)= locpoles(pp);
 210     }
 211
 212     locqq = 0.0;
 213     locdqq = 0.0;
 214     for (pp = 1; pp <= deg; pp++) {
 215       Inverse =  1.0 / (flatknots(kindex + pp)  - flatknots(kindex - qq + pp + 1));
 216       Saved = (U - flatknots(kindex - qq + pp + 1)) * Inverse * locpoles(pp);
 217       locpoles(pp) *= (flatknots(kindex + pp) - U) * Inverse;
 218       locpoles(pp) += locqq;
 219       locqq = Saved ;
 220       LocalInverse = (Standard_Real) (deg) * Inverse;
 221       Saved = LocalInverse * locdpoles(pp);
 222       locdpoles(pp) *= - LocalInverse ;
 223       locdpoles(pp) +=   locdqq;
 224       locdqq = Saved ;
 225     }
 226
 227     locpoles(qq) = locqq;
 228     locdpoles(qq) = locdqq;
 229
 230     for (j = 1; j <= deg1; j++) {
 231       val = locpoles(j);
 232       theindex = j+oldkindex-deg1;
 233       A(i, theindex) = val;
 234       DA(i, theindex) = locdpoles(j);
 235     }
 236
 237     for (j = 1; j < oldkindex-deg; j++)
 238       A(i, j) = DA(i, j) = 0.0;
 239     for (j = oldkindex+1; j <= nbpoles; j++)
 240       A(i, j) = DA(i, j) = 0.0;
 241
 242   }
 243
 244 }