src/Law/Law.cxx

   1 // Created on: 1995-11-20
   2 // Created by: Laurent BOURESCHE
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17
  18 #include <Adaptor3d_Curve.hxx>
  19 #include <BSplCLib.hxx>
  20 #include <gp_Pnt.hxx>
  21 #include <Law.hxx>
  22 #include <Law_BSpFunc.hxx>
  23 #include <Law_BSpline.hxx>
  24 #include <Law_Interpolate.hxx>
  25 #include <Law_Linear.hxx>
  26 #include <TColStd_Array1OfBoolean.hxx>
  27 #include <TColStd_Array1OfReal.hxx>
  28 #include <TColStd_HArray1OfBoolean.hxx>
  29 #include <TColStd_HArray1OfReal.hxx>
  30
  31 Handle(Law_BSpFunc) Law::MixBnd(const Handle(Law_Linear)& Lin)
  32 {
  33   Standard_Real f,l;
  34   Lin->Bounds(f,l);
  35   TColStd_Array1OfReal Knots(1,4);
  36   TColStd_Array1OfInteger Mults(1,4);
  37   Knots(1) = f; Knots(4) = l;
  38   Knots(2) = 0.75*f + 0.25*l; Knots(3) = 0.25*f + 0.75*l;
  39   Mults(1) = Mults(4) = 4;
  40   Mults(2) = Mults(3) = 1;
  41   Handle(TColStd_HArray1OfReal) pol = Law::MixBnd(3,Knots,Mults,Lin);
  42   Handle(Law_BSpline) bs = new Law_BSpline(pol->Array1(),Knots,Mults,3);
  43   Handle(Law_BSpFunc) bsf = new Law_BSpFunc();
  44   bsf->SetCurve(bs);
  45   return bsf;
  46 }
  47
  48 Handle(TColStd_HArray1OfReal) Law::MixBnd
  49 (const Standard_Integer         Degree,
  50  const TColStd_Array1OfReal&    Knots,
  51  const TColStd_Array1OfInteger& Mults,
  52  const Handle(Law_Linear)&      Lin)
  53 {
  54   Standard_Integer nbpol = 0, nbfk = 0, i, j, k = 0;
  55   for(i = Mults.Lower(); i <= Mults.Upper(); i++){
  56     nbfk += Mults(i);
  57   }
  58   TColStd_Array1OfReal fk(1,nbfk);
  59   for(i = Mults.Lower(); i <= Mults.Upper(); i++){
  60     for(j = 1; j <= Mults(i); j++){
  61       fk(++k) = Knots(i);
  62     }
  63   }
  64   nbpol = nbfk - Degree - 1;
  65   TColStd_Array1OfReal par(1,nbpol);
  66   BSplCLib::BuildSchoenbergPoints(Degree,fk,par);
  67   Handle(TColStd_HArray1OfReal) res = new TColStd_HArray1OfReal(1,nbpol);
  68   TColStd_Array1OfReal& pol = res->ChangeArray1();
  69   for(i = 1; i <= nbpol; i++){
  70     pol(i) = Lin->Value(par(i));
  71   }
  72   TColStd_Array1OfInteger ord(1,nbpol);
  73   ord.Init(0);
  74   BSplCLib::Interpolate(Degree,fk,par,ord,1,pol(1),i);
  75   if(nbpol >= 4){
  76     pol(2) = pol(1);
  77     pol(nbpol - 1) = pol(nbpol);
  78   }
  79   return res;
  80 }
  81
  82 static Standard_Real eval1(const Standard_Real    p,
  83                            const Standard_Real    first,
  84                            const Standard_Real    last,
  85                            const Standard_Real    piv,
  86                            const Standard_Boolean nulr)
  87 {
  88   if((nulr && p >= piv) || (!nulr && p <= piv)) return 0.;
  89   else if(nulr){
  90     Standard_Real a = piv - first;
  91     a *= a;
  92     a = 1./a;
  93     Standard_Real b = p - first;
  94     a *= b;
  95     b = piv - p;
  96     a *= b;
  97     a *= b;
  98     return a;
  99   }
 100   else{
 101     Standard_Real a = last - piv;
 102     a *= a;
 103     a = 1./a;
 104     Standard_Real b = last - p;
 105     a *= b;
 106     b = p - piv;
 107     a *= b;
 108     a *= b;
 109     return a;
 110   }
 111 }
 112
 113 Handle(TColStd_HArray1OfReal) Law::MixTgt
 114 (const Standard_Integer         Degree,
 115  const TColStd_Array1OfReal&    Knots,
 116  const TColStd_Array1OfInteger& Mults,
 117  const Standard_Boolean         NulOnTheRight,
 118  const Standard_Integer         Index)
 119 {
 120   Standard_Real first = Knots(Knots.Lower());
 121   Standard_Real last = Knots(Knots.Upper());
 122   Standard_Real piv = Knots(Index);
 123   Standard_Integer nbpol = 0, nbfk = 0, i, j, k = 0;
 124   for(i = Mults.Lower(); i <= Mults.Upper(); i++){
 125     nbfk += Mults(i);
 126   }
 127   TColStd_Array1OfReal fk(1,nbfk);
 128   for(i = Mults.Lower(); i <= Mults.Upper(); i++){
 129     for(j = 1; j <= Mults(i); j++){
 130       fk(++k) = Knots(i);
 131     }
 132   }
 133   nbpol = nbfk - Degree - 1;
 134   TColStd_Array1OfReal par(1,nbpol);
 135   BSplCLib::BuildSchoenbergPoints(Degree,fk,par);
 136   Handle(TColStd_HArray1OfReal) res = new TColStd_HArray1OfReal(1,nbpol);
 137   TColStd_Array1OfReal& pol = res->ChangeArray1();
 138   for(i = 1; i <= nbpol; i++){
 139     pol(i) = eval1(par(i),first,last,piv,NulOnTheRight);
 140   }
 141   TColStd_Array1OfInteger ord(1,nbpol);
 142   ord.Init(0);
 143   BSplCLib::Interpolate(Degree,fk,par,ord,1,pol(1),i);
 144   return res;
 145 }
 146
 147 Handle(Law_BSpline) Law::Reparametrize(const Adaptor3d_Curve&   Curve,
 148                                        const Standard_Real    First,
 149                                        const Standard_Real    Last,
 150                                        const Standard_Boolean HasDF,
 151                                        const Standard_Boolean HasDL,
 152                                        const Standard_Real    DFirst,
 153                                        const Standard_Real    DLast,
 154                                        const Standard_Boolean Rev,
 155                                        const Standard_Integer NbPoints)
 156 {
 157   // On evalue la longeur approximative de la courbe.
 158
 159   Standard_Integer i;
 160   Standard_Real DDFirst = DFirst, DDLast = DLast;
 161   if(HasDF && Rev) DDFirst = -DFirst;
 162   if(HasDL && Rev) DDLast = -DLast;
 163   TColStd_Array1OfReal cumdist(1,2*NbPoints);
 164   TColStd_Array1OfReal ucourbe(1,2*NbPoints);
 165   gp_Pnt P1,P2;
 166   Standard_Real U1 = Curve.FirstParameter();
 167   Standard_Real U2 = Curve.LastParameter();
 168   Standard_Real U, DU, Length = 0.;
 169   if(!Rev){
 170     P1 = Curve.Value(U1);
 171     U = U1;
 172     DU = (U2 - U1) / ( 2*NbPoints - 1);
 173   }
 174   else{
 175     P1 = Curve.Value(U2);
 176     U = U2;
 177     DU = (U1 - U2) / ( 2*NbPoints - 1);
 178   }
 179   for ( i = 1; i <= 2*NbPoints ; i ++) {
 180     P2 = Curve.Value(U);
 181     Length += P2.Distance(P1);
 182     cumdist(i) = Length;
 183     ucourbe(i) = U;
 184     U += DU;
 185     P1 = P2;
 186   }
 187   if(Rev) ucourbe(2*NbPoints) = U1;
 188   else ucourbe(2*NbPoints) = U2;
 189
 190   Handle(TColStd_HArray1OfReal) point = new TColStd_HArray1OfReal(1,NbPoints);
 191   Handle(TColStd_HArray1OfReal) param = new TColStd_HArray1OfReal(1,NbPoints);
 192
 193   Standard_Real DCorde = Length / ( NbPoints - 1);
 194   Standard_Real Corde  = DCorde;
 195   Standard_Integer Index = 1;
 196   Standard_Real Alpha;
 197   Standard_Real fac = 1./(NbPoints-1);
 198
 199   point->SetValue(1,ucourbe(1));
 200   param->SetValue(1,First);
 201   point->SetValue(NbPoints,ucourbe(2*NbPoints));
 202   param->SetValue(NbPoints,Last);
 203
 204   for ( i = 2; i < NbPoints; i++) {
 205     while (cumdist(Index) < Corde) Index ++;
 206
 207     Alpha = (Corde - cumdist(Index - 1)) / (cumdist(Index) -cumdist(Index - 1));
 208     U = ucourbe(Index-1) + Alpha * (ucourbe(Index) - ucourbe(Index-1));
 209     point->SetValue(i,U);
 210     param->SetValue(i,((NbPoints-i)*First+(i-1)*Last)*fac);
 211     Corde = i*DCorde;
 212   }
 213   Law_Interpolate inter(point,param,0,1.e-9);
 214   if(HasDF || HasDL){
 215     TColStd_Array1OfReal tgt(1,NbPoints);
 216     Handle(TColStd_HArray1OfBoolean)
 217       flag = new TColStd_HArray1OfBoolean(1,NbPoints);
 218     flag->ChangeArray1().Init(0);
 219     if(HasDF){ flag->SetValue(1,1); tgt.SetValue(1,DDFirst); }
 220     if(HasDL){ flag->SetValue(NbPoints,1); tgt.SetValue(NbPoints,DDLast); }
 221     inter.Load(tgt,flag);
 222   }
 223   inter.Perform();
 224   if(!inter.IsDone())
 225     throw Standard_Failure("Law::Reparametrize echec interpolation");
 226   Handle(Law_BSpline) bs = inter.Curve();
 227   return bs;
 228 }
 229
 230 static Standard_Real eval2(const Standard_Real        p,
 231 //                         const Standard_Real        first,
 232                            const Standard_Real        ,
 233 //                         const Standard_Real        last,
 234                            const Standard_Real        ,
 235                            const Standard_Real        mil,
 236                            const Standard_Boolean     hasfirst,
 237                            const Standard_Boolean     haslast,
 238                            const Handle(Law_BSpline)& bs1,
 239                            const Handle(Law_BSpline)& bs2)
 240 {
 241   if(hasfirst && p < mil) return bs1->Value(p);
 242   else if(haslast && p > mil) return bs2->Value(p);
 243   else return 1.;
 244 }
 245
 246 Handle(Law_BSpline) Law::Scale(const Standard_Real    First,
 247                                const Standard_Real    Last,
 248                                const Standard_Boolean HasF,
 249                                const Standard_Boolean HasL,
 250                                const Standard_Real    VFirst,
 251                                const Standard_Real    VLast)
 252 {
 253   Standard_Integer i;
 254   Standard_Real Milieu = 0.5*(First+Last);
 255   TColStd_Array1OfReal   knot(1,3);
 256   TColStd_Array1OfReal   fknot(1,10);
 257   TColStd_Array1OfInteger mult(1,3);
 258   knot(1) = First; knot(2) = Milieu; knot(3) = Last;
 259   fknot(1) = fknot(2) = fknot(3) = fknot(4) = First;
 260   fknot(10) = fknot(9) = fknot(8) = fknot(7) = Last;
 261   fknot(5) = fknot(6) =  Milieu;
 262   mult(1) = 4; mult(3) = 4; mult(2) = 2;
 263
 264   TColStd_Array1OfReal pbs(1,4);
 265   TColStd_Array1OfReal kbs(1,2);
 266   TColStd_Array1OfInteger mbs(1,2);
 267   mbs(1) = mbs(2) = 4;
 268   Handle(Law_BSpline) bs1,bs2;
 269   if(HasF){
 270     pbs(1) = pbs(2) = VFirst;
 271     pbs(3) = pbs(4) = 1.;
 272     kbs(1) = First; kbs(2) = Milieu;
 273     bs1 = new Law_BSpline(pbs,kbs,mbs,3);
 274   }
 275   if(HasL){
 276     pbs(1) = pbs(2) = 1.;
 277     pbs(3) = pbs(4) = VLast;
 278     kbs(1) = Milieu; kbs(2) = Last;
 279     bs2 = new Law_BSpline(pbs,kbs,mbs,3);
 280   }
 281
 282   TColStd_Array1OfReal pol(1,6);
 283   TColStd_Array1OfReal par(1,6);
 284   BSplCLib::BuildSchoenbergPoints(3,fknot,par);
 285   for(i = 1; i <= 6; i++){
 286     pol(i) = eval2(par(i),First,Last,Milieu,HasF,HasL,bs1,bs2);
 287   }
 288   TColStd_Array1OfInteger ord(1,6);
 289   ord.Init(0);
 290   BSplCLib::Interpolate(3,fknot,par,ord,1,pol(1),i);
 291   bs1 = new Law_BSpline(pol,knot,mult,3);
 292   return bs1;
 293 }
 294
 295
 296 Handle(Law_BSpline) Law::ScaleCub(const Standard_Real    First,
 297                                   const Standard_Real    Last,
 298                                   const Standard_Boolean HasF,
 299                                   const Standard_Boolean HasL,
 300                                   const Standard_Real    VFirst,
 301                                   const Standard_Real    VLast)
 302 {
 303 //Standard_Integer i;
 304   Standard_Real Milieu = 0.5*(First+Last);
 305   TColStd_Array1OfReal pol(1,5);
 306   TColStd_Array1OfReal   knot(1,3);
 307   TColStd_Array1OfInteger mult(1,3);
 308   knot(1) = First; knot(2) = Milieu; knot(3) = Last;
 309   mult(1) = 4; mult(3) = 4; mult(2) = 1;
 310
 311   Handle(Law_BSpline) bs;
 312   if(HasF){ pol(1) = pol(2) = VFirst;}
 313   else { pol(1) = pol(2) = 1.; }
 314   if(HasL){ pol(4) = pol(5) = VLast;}
 315   else { pol(4) = pol(5) = 1.; }
 316
 317   pol(3) = 1.;
 318
 319   bs = new Law_BSpline(pol,knot,mult,3);
 320   return bs;
 321 }