| 1 | // Created on: 1993-07-07 |
| 2 | // Created by: Jean Claude VAUTHIER |
| 3 | // Copyright (c) 1993-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 | // Version: |
| 18 | //pmn 24/09/96 Ajout du prolongement de courbe. |
| 19 | // jct 15/04/97 Ajout du prolongement de surface. |
| 20 | // jct 24/04/97 simplification ou suppression de calculs |
| 21 | // inutiles dans ExtendSurfByLength |
| 22 | // correction de Tbord et Continuity=0 accepte |
| 23 | // correction du calcul de lambda et appel a |
| 24 | // TangExtendToConstraint avec lambmin au lieu de 1. |
| 25 | // correction du passage Sr rat --> BSp nD |
| 26 | // xab 26/06/97 treatement partiel anulation des derivees |
| 27 | // partiels du denonimateur des Surfaces BSplines Rationnelles |
| 28 | // dans le cas de valeurs proportionnelles des denominateurs |
| 29 | // en umin umax et/ou vmin vmax. |
| 30 | // pmn 4/07/97 Gestion de la continuite dans BuildCurve3d (PRO9097) |
| 31 | |
| 32 | // xab 10/07/97 on revient en arriere sur l'ajout du 26/06/97 |
| 33 | // pmn 26/09/97 Ajout des parametres d'approx dans BuildCurve3d |
| 34 | // xab 29/09/97 on reintegre l'ajout du 26/06/97 |
| 35 | // pmn 31/10/97 Ajoute AdjustExtremity |
| 36 | // jct 26/11/98 blindage dans ExtendSurf qd NTgte = 0 (CTS21288) |
| 37 | // jct 19/01/99 traitement de la periodicite dans ExtendSurf |
| 38 | // Design: |
| 39 | // Warning: None |
| 40 | // References: None |
| 41 | // Language: C++2.0 |
| 42 | // Purpose: |
| 43 | |
| 44 | // Declarations: |
| 45 | |
| 46 | #include <GeomLib.ixx> |
| 47 | |
| 48 | #include <Precision.hxx> |
| 49 | #include <GeomConvert.hxx> |
| 50 | #include <Hermit.hxx> |
| 51 | #include <Standard_NotImplemented.hxx> |
| 52 | #include <GeomLib_MakeCurvefromApprox.hxx> |
| 53 | #include <GeomLib_DenominatorMultiplier.hxx> |
| 54 | #include <GeomLib_DenominatorMultiplierPtr.hxx> |
| 55 | #include <GeomLib_PolyFunc.hxx> |
| 56 | #include <GeomLib_LogSample.hxx> |
| 57 | |
| 58 | #include <AdvApprox_ApproxAFunction.hxx> |
| 59 | #include <AdvApprox_PrefAndRec.hxx> |
| 60 | |
| 61 | #include <Adaptor2d_HCurve2d.hxx> |
| 62 | #include <Adaptor3d_HCurve.hxx> |
| 63 | #include <Adaptor3d_HSurface.hxx> |
| 64 | #include <Adaptor3d_CurveOnSurface.hxx> |
| 65 | #include <Geom2dAdaptor_Curve.hxx> |
| 66 | #include <GeomAdaptor_Surface.hxx> |
| 67 | #include <GeomAdaptor_HSurface.hxx> |
| 68 | #include <Geom2dAdaptor_HCurve.hxx> |
| 69 | #include <Geom2dAdaptor_GHCurve.hxx> |
| 70 | |
| 71 | #include <Geom2d_BSplineCurve.hxx> |
| 72 | #include <Geom_BSplineCurve.hxx> |
| 73 | #include <Geom2d_BezierCurve.hxx> |
| 74 | #include <Geom_BezierCurve.hxx> |
| 75 | #include <Geom_RectangularTrimmedSurface.hxx> |
| 76 | #include <Geom_Plane.hxx> |
| 77 | #include <Geom_Line.hxx> |
| 78 | #include <Geom2d_Line.hxx> |
| 79 | #include <Geom_Circle.hxx> |
| 80 | #include <Geom2d_Circle.hxx> |
| 81 | #include <Geom_Ellipse.hxx> |
| 82 | #include <Geom2d_Ellipse.hxx> |
| 83 | #include <Geom_Parabola.hxx> |
| 84 | #include <Geom2d_Parabola.hxx> |
| 85 | #include <Geom_Hyperbola.hxx> |
| 86 | #include <Geom2d_Hyperbola.hxx> |
| 87 | #include <Geom_TrimmedCurve.hxx> |
| 88 | #include <Geom2d_TrimmedCurve.hxx> |
| 89 | #include <Geom_OffsetCurve.hxx> |
| 90 | #include <Geom2d_OffsetCurve.hxx> |
| 91 | #include <Geom_BezierSurface.hxx> |
| 92 | #include <Geom_BSplineSurface.hxx> |
| 93 | |
| 94 | #include <BSplCLib.hxx> |
| 95 | #include <BSplSLib.hxx> |
| 96 | #include <PLib.hxx> |
| 97 | #include <math_Matrix.hxx> |
| 98 | #include <math_Vector.hxx> |
| 99 | #include <math_Jacobi.hxx> |
| 100 | #include <math.hxx> |
| 101 | #include <math_FunctionAllRoots.hxx> |
| 102 | #include <math_FunctionSample.hxx> |
| 103 | |
| 104 | #include <TColStd_HArray1OfReal.hxx> |
| 105 | #include <TColgp_Array1OfPnt.hxx> |
| 106 | #include <TColgp_Array1OfVec.hxx> |
| 107 | #include <TColgp_Array2OfPnt.hxx> |
| 108 | #include <TColgp_HArray2OfPnt.hxx> |
| 109 | #include <TColgp_Array1OfPnt2d.hxx> |
| 110 | #include <TColgp_Array1OfXYZ.hxx> |
| 111 | #include <TColStd_Array1OfReal.hxx> |
| 112 | #include <TColStd_Array2OfReal.hxx> |
| 113 | #include <TColStd_HArray2OfReal.hxx> |
| 114 | #include <TColStd_Array1OfInteger.hxx> |
| 115 | |
| 116 | #include <gp_TrsfForm.hxx> |
| 117 | #include <gp_Lin.hxx> |
| 118 | #include <gp_Lin2d.hxx> |
| 119 | #include <gp_Circ.hxx> |
| 120 | #include <gp_Circ2d.hxx> |
| 121 | #include <gp_Elips.hxx> |
| 122 | #include <gp_Elips2d.hxx> |
| 123 | #include <gp_Hypr.hxx> |
| 124 | #include <gp_Hypr2d.hxx> |
| 125 | #include <gp_Parab.hxx> |
| 126 | #include <gp_Parab2d.hxx> |
| 127 | #include <gp_GTrsf2d.hxx> |
| 128 | #include <gp_Trsf2d.hxx> |
| 129 | |
| 130 | #include <ElCLib.hxx> |
| 131 | #include <Geom2dConvert.hxx> |
| 132 | #include <GeomConvert_CompCurveToBSplineCurve.hxx> |
| 133 | #include <GeomConvert_ApproxSurface.hxx> |
| 134 | |
| 135 | #include <CSLib.hxx> |
| 136 | #include <CSLib_NormalStatus.hxx> |
| 137 | |
| 138 | |
| 139 | #include <Standard_ConstructionError.hxx> |
| 140 | |
| 141 | //======================================================================= |
| 142 | //function : ComputeLambda |
| 143 | //purpose : Calcul le facteur lambda qui minimise la variation de vittesse |
| 144 | // sur une interpolation d'hermite d'ordre (i,0) |
| 145 | //======================================================================= |
| 146 | static void ComputeLambda(const math_Matrix& Constraint, |
| 147 | const math_Matrix& Hermit, |
| 148 | const Standard_Real Length, |
| 149 | Standard_Real& Lambda ) |
| 150 | { |
| 151 | Standard_Integer size = Hermit.RowNumber(); |
| 152 | Standard_Integer Continuity = size-2; |
| 153 | Standard_Integer ii, jj, ip, pp; |
| 154 | |
| 155 | //Minimization |
| 156 | math_Matrix HDer(1, size-1, 1, size); |
| 157 | for (jj=1; jj<=size; jj++) { |
| 158 | for (ii=1; ii<size;ii++) { |
| 159 | HDer(ii, jj) = ii*Hermit(jj, ii+1); |
| 160 | } |
| 161 | } |
| 162 | |
| 163 | math_Vector V(1, size); |
| 164 | math_Vector Vec1(1, Constraint.RowNumber()); |
| 165 | math_Vector Vec2(1, Constraint.RowNumber()); |
| 166 | math_Vector Vec3(1, Constraint.RowNumber()); |
| 167 | math_Vector Vec4(1, Constraint.RowNumber()); |
| 168 | |
| 169 | Standard_Real * polynome = &HDer(1,1); |
| 170 | Standard_Real * valhder = &V(1); |
| 171 | Vec2 = Constraint.Col(2); |
| 172 | Vec2 /= Length; |
| 173 | Standard_Real t, squared1 = Vec2.Norm2(), GW; |
| 174 | // math_Matrix Vec(1, Constraint.RowNumber(), 1, size-1); |
| 175 | // gp_Vec Vfirst(p0.XYZ()), Vlast(Point.XYZ()); |
| 176 | // TColgp_Array1OfVec Der(2, 4); |
| 177 | // Der(2) = d1; Der(3) = d2; Der(4) = d3; |
| 178 | |
| 179 | Standard_Integer GOrdre = 4 + 4*Continuity, |
| 180 | DDim=Continuity*(Continuity+2); |
| 181 | math_Vector GaussP(1, GOrdre), GaussW(1, GOrdre), |
| 182 | pol2(1, 2*Continuity+1), |
| 183 | pol4(1, 4*Continuity+1); |
| 184 | math::GaussPoints(GOrdre, GaussP); |
| 185 | math::GaussWeights (GOrdre, GaussW); |
| 186 | pol4.Init(0.); |
| 187 | |
| 188 | for (ip=1; ip<=GOrdre; ip++) { |
| 189 | t = (GaussP(ip)+1.)/2; |
| 190 | GW = GaussW(ip); |
| 191 | PLib::NoDerivativeEvalPolynomial(t , Continuity, Continuity+2, DDim, |
| 192 | polynome[0], valhder[0]); |
| 193 | V /= Length; //Normalisation |
| 194 | |
| 195 | // i |
| 196 | // C'(t) = SUM Vi*Lambda |
| 197 | Vec1 = Constraint.Col(1); |
| 198 | Vec1 *= V(1); |
| 199 | Vec1 += V(size)*Constraint.Col(size); |
| 200 | Vec2 = Constraint.Col(2); |
| 201 | Vec2 *= V(2); |
| 202 | if (Continuity > 1) { |
| 203 | Vec3 = Constraint.Col(3); |
| 204 | Vec3 *= V(3); |
| 205 | if (Continuity > 2) { |
| 206 | Vec4 = Constraint.Col(4); |
| 207 | Vec4 *= V(4); |
| 208 | } |
| 209 | } |
| 210 | |
| 211 | |
| 212 | // 2 2 |
| 213 | // C'(t) - C'(0) |
| 214 | |
| 215 | pol2(1) = Vec1.Norm2(); |
| 216 | pol2(2) = 2*(Vec1.Multiplied(Vec2)); |
| 217 | pol2(3) = Vec2.Norm2() - squared1; |
| 218 | if (Continuity>1) { |
| 219 | pol2(3) += 2*(Vec1.Multiplied(Vec3)); |
| 220 | pol2(4) = 2*(Vec2.Multiplied(Vec3)); |
| 221 | pol2(5) = Vec3.Norm2(); |
| 222 | if (Continuity>2) { |
| 223 | pol2(4)+= 2*(Vec1.Multiplied(Vec4)); |
| 224 | pol2(5)+= 2*(Vec2.Multiplied(Vec4)); |
| 225 | pol2(6) = 2*(Vec3.Multiplied(Vec4)); |
| 226 | pol2(7) = Vec4.Norm2(); |
| 227 | } |
| 228 | } |
| 229 | |
| 230 | // 2 2 2 |
| 231 | // Integrale de ( C'(t) - C'(0) ) |
| 232 | for (ii=1; ii<=pol2.Length(); ii++) { |
| 233 | pp = ii; |
| 234 | for(jj=1; jj<ii; jj++, pp++) { |
| 235 | pol4(pp) += 2*GW*pol2(ii)*pol2(jj); |
| 236 | } |
| 237 | pol4(2*ii-1) += GW*Pow(pol2(ii), 2); |
| 238 | } |
| 239 | } |
| 240 | |
| 241 | Standard_Real EMin, E; |
| 242 | PLib::NoDerivativeEvalPolynomial(Lambda , pol4.Length()-1, 1, |
| 243 | pol4.Length()-1, |
| 244 | pol4(1), EMin); |
| 245 | |
| 246 | if (EMin > Precision::Confusion()) { |
| 247 | // Recheche des extrema de la fonction |
| 248 | GeomLib_PolyFunc FF(pol4); |
| 249 | GeomLib_LogSample S(Lambda/1000, 50*Lambda, 100); |
| 250 | math_FunctionAllRoots Solve(FF, S, Precision::Confusion(), |
| 251 | Precision::Confusion()*(Length+1), |
| 252 | 1.e-15); |
| 253 | if (Solve.IsDone()) { |
| 254 | for (ii=1; ii<=Solve.NbPoints(); ii++) { |
| 255 | t = Solve.GetPoint(ii); |
| 256 | PLib::NoDerivativeEvalPolynomial(t , pol4.Length()-1, 1, |
| 257 | pol4.Length()-1, |
| 258 | pol4(1), E); |
| 259 | if (E < EMin) { |
| 260 | Lambda = t; |
| 261 | EMin = E; |
| 262 | } |
| 263 | } |
| 264 | } |
| 265 | } |
| 266 | } |
| 267 | |
| 268 | #include <Extrema_LocateExtPC.hxx> |
| 269 | //======================================================================= |
| 270 | //function : RemovePointsFromArray |
| 271 | //purpose : |
| 272 | //======================================================================= |
| 273 | |
| 274 | void GeomLib::RemovePointsFromArray(const Standard_Integer NumPoints, |
| 275 | const TColStd_Array1OfReal& InParameters, |
| 276 | Handle(TColStd_HArray1OfReal)& OutParameters) |
| 277 | { |
| 278 | Standard_Integer ii, |
| 279 | jj, |
| 280 | add_one_point, |
| 281 | loc_num_points, |
| 282 | num_points, |
| 283 | index ; |
| 284 | Standard_Real delta, |
| 285 | current_parameter ; |
| 286 | |
| 287 | loc_num_points = Max(0,NumPoints-2) ; |
| 288 | delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ; |
| 289 | delta /= (Standard_Real) (loc_num_points + 1) ; |
| 290 | num_points = 1 ; |
| 291 | current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ; |
| 292 | ii = InParameters.Lower() + 1 ; |
| 293 | for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) { |
| 294 | add_one_point = 0 ; |
| 295 | while ( ii < InParameters.Upper() && InParameters(ii) < current_parameter) { |
| 296 | ii += 1 ; |
| 297 | add_one_point = 1 ; |
| 298 | } |
| 299 | num_points += add_one_point ; |
| 300 | current_parameter += delta ; |
| 301 | } |
| 302 | if (NumPoints <= 2) { |
| 303 | num_points = 2 ; |
| 304 | } |
| 305 | index = 2 ; |
| 306 | current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ; |
| 307 | OutParameters = |
| 308 | new TColStd_HArray1OfReal(1,num_points) ; |
| 309 | OutParameters->ChangeArray1()(1) = InParameters(InParameters.Lower()) ; |
| 310 | ii = InParameters.Lower() + 1 ; |
| 311 | for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) { |
| 312 | add_one_point = 0 ; |
| 313 | while (ii < InParameters.Upper() && InParameters(ii) < current_parameter) { |
| 314 | ii += 1 ; |
| 315 | add_one_point = 1 ; |
| 316 | } |
| 317 | if (add_one_point && index <= num_points) { |
| 318 | OutParameters->ChangeArray1()(index) = InParameters(ii-1) ; |
| 319 | index += 1 ; |
| 320 | } |
| 321 | current_parameter += delta ; |
| 322 | } |
| 323 | OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ; |
| 324 | } |
| 325 | //======================================================================= |
| 326 | //function : DensifyArray1OfReal |
| 327 | //purpose : |
| 328 | //======================================================================= |
| 329 | |
| 330 | void GeomLib::DensifyArray1OfReal(const Standard_Integer MinNumPoints, |
| 331 | const TColStd_Array1OfReal& InParameters, |
| 332 | Handle(TColStd_HArray1OfReal)& OutParameters) |
| 333 | { |
| 334 | Standard_Integer ii, |
| 335 | in_order, |
| 336 | num_points, |
| 337 | num_parameters_to_add, |
| 338 | index ; |
| 339 | Standard_Real delta, |
| 340 | current_parameter ; |
| 341 | |
| 342 | in_order = 1 ; |
| 343 | if (MinNumPoints > InParameters.Length()) { |
| 344 | |
| 345 | // |
| 346 | // checks the paramaters are in increasing order |
| 347 | // |
| 348 | for (ii = InParameters.Lower() ; ii < InParameters.Upper() ; ii++) { |
| 349 | if (InParameters(ii) > InParameters(ii+1)) { |
| 350 | in_order = 0 ; |
| 351 | break ; |
| 352 | } |
| 353 | } |
| 354 | if (in_order) { |
| 355 | num_parameters_to_add = MinNumPoints - InParameters.Length() ; |
| 356 | delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ; |
| 357 | delta /= (Standard_Real) (num_parameters_to_add + 1) ; |
| 358 | num_points = MinNumPoints ; |
| 359 | OutParameters = |
| 360 | new TColStd_HArray1OfReal(1,num_points) ; |
| 361 | index = 1 ; |
| 362 | current_parameter = InParameters(InParameters.Lower()) ; |
| 363 | OutParameters->ChangeArray1()(index) = current_parameter ; |
| 364 | index += 1 ; |
| 365 | current_parameter += delta ; |
| 366 | for (ii = InParameters.Lower() + 1 ; index <= num_points && ii <= InParameters.Upper() ; ii++) { |
| 367 | while (current_parameter < InParameters(ii) && index <= num_points) { |
| 368 | OutParameters->ChangeArray1()(index) = current_parameter ; |
| 369 | index += 1 ; |
| 370 | current_parameter += delta ; |
| 371 | } |
| 372 | if (index <= num_points) { |
| 373 | OutParameters->ChangeArray1()(index) = InParameters(ii) ; |
| 374 | } |
| 375 | index += 1 ; |
| 376 | } |
| 377 | // |
| 378 | // beware of roundoff ! |
| 379 | // |
| 380 | OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ; |
| 381 | } |
| 382 | else { |
| 383 | index = 1 ; |
| 384 | num_points = InParameters.Length() ; |
| 385 | OutParameters = |
| 386 | new TColStd_HArray1OfReal(1,num_points) ; |
| 387 | for (ii = InParameters.Lower() ; ii <= InParameters.Upper() ; ii++) { |
| 388 | OutParameters->ChangeArray1()(index) = InParameters(ii) ; |
| 389 | index += 1 ; |
| 390 | } |
| 391 | } |
| 392 | } |
| 393 | else { |
| 394 | index = 1 ; |
| 395 | num_points = InParameters.Length() ; |
| 396 | OutParameters = |
| 397 | new TColStd_HArray1OfReal(1,num_points) ; |
| 398 | for (ii = InParameters.Lower() ; ii <= InParameters.Upper() ; ii++) { |
| 399 | OutParameters->ChangeArray1()(index) = InParameters(ii) ; |
| 400 | index += 1 ; |
| 401 | } |
| 402 | } |
| 403 | } |
| 404 | |
| 405 | //======================================================================= |
| 406 | //function : FuseIntervals |
| 407 | //purpose : |
| 408 | //======================================================================= |
| 409 | void GeomLib::FuseIntervals(const TColStd_Array1OfReal& I1, |
| 410 | const TColStd_Array1OfReal& I2, |
| 411 | TColStd_SequenceOfReal& Seq, |
| 412 | const Standard_Real Epspar) |
| 413 | { |
| 414 | Standard_Integer ind1=1, ind2=1; |
| 415 | Standard_Real v1, v2; |
| 416 | // Initialisations : les IND1 et IND2 pointent sur le 1er element |
| 417 | // de chacune des 2 tables a traiter.INDS pointe sur le dernier |
| 418 | // element cree de TABSOR |
| 419 | |
| 420 | |
| 421 | //--- On remplit TABSOR en parcourant TABLE1 et TABLE2 simultanement --- |
| 422 | //------------------ en eliminant les occurrences multiples ------------ |
| 423 | |
| 424 | while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) { |
| 425 | v1 = I1(ind1); |
| 426 | v2 = I2(ind2); |
| 427 | if (Abs(v1-v2)<= Epspar) { |
| 428 | // Ici les elements de I1 et I2 conviennent . |
| 429 | Seq.Append((v1+v2)/2); |
| 430 | ind1++; |
| 431 | ind2++; |
| 432 | } |
| 433 | else if (v1 < v2) { |
| 434 | // Ici l' element de I1 convient. |
| 435 | Seq.Append(v1); |
| 436 | ind1++; |
| 437 | } |
| 438 | else { |
| 439 | // Ici l' element de TABLE2 convient. |
| 440 | Seq.Append(v2); |
| 441 | ind2++; |
| 442 | } |
| 443 | } |
| 444 | |
| 445 | if (ind1>I1.Upper()) { |
| 446 | //----- Ici I1 est epuise, on complete avec la fin de TABLE2 ------- |
| 447 | |
| 448 | for (; ind2<=I2.Upper(); ind2++) { |
| 449 | Seq.Append(I2(ind2)); |
| 450 | } |
| 451 | } |
| 452 | |
| 453 | if (ind2>I2.Upper()) { |
| 454 | //----- Ici I2 est epuise, on complete avec la fin de I1 ------- |
| 455 | for (; ind1<=I1.Upper(); ind1++) { |
| 456 | Seq.Append(I1(ind1)); |
| 457 | } |
| 458 | } |
| 459 | } |
| 460 | |
| 461 | |
| 462 | //======================================================================= |
| 463 | //function : EvalMaxParametricDistance |
| 464 | //purpose : |
| 465 | //======================================================================= |
| 466 | |
| 467 | void GeomLib::EvalMaxParametricDistance(const Adaptor3d_Curve& ACurve, |
| 468 | const Adaptor3d_Curve& AReferenceCurve, |
| 469 | // const Standard_Real Tolerance, |
| 470 | const Standard_Real , |
| 471 | const TColStd_Array1OfReal& Parameters, |
| 472 | Standard_Real& MaxDistance) |
| 473 | { |
| 474 | Standard_Integer ii ; |
| 475 | |
| 476 | Standard_Real max_squared = 0.0e0, |
| 477 | // tolerance_squared, |
| 478 | local_distance_squared ; |
| 479 | |
| 480 | // tolerance_squared = Tolerance * Tolerance ; |
| 481 | gp_Pnt Point1 ; |
| 482 | gp_Pnt Point2 ; |
| 483 | for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) { |
| 484 | ACurve.D0(Parameters(ii), |
| 485 | Point1) ; |
| 486 | AReferenceCurve.D0(Parameters(ii), |
| 487 | Point2) ; |
| 488 | local_distance_squared = |
| 489 | Point1.SquareDistance (Point2) ; |
| 490 | max_squared = Max(max_squared,local_distance_squared) ; |
| 491 | } |
| 492 | if (max_squared > 0.0e0) { |
| 493 | MaxDistance = sqrt(max_squared) ; |
| 494 | } |
| 495 | else { |
| 496 | MaxDistance = 0.0e0 ; |
| 497 | } |
| 498 | |
| 499 | } |
| 500 | //======================================================================= |
| 501 | //function : EvalMaxDistanceAlongParameter |
| 502 | //purpose : |
| 503 | //======================================================================= |
| 504 | |
| 505 | void GeomLib::EvalMaxDistanceAlongParameter(const Adaptor3d_Curve& ACurve, |
| 506 | const Adaptor3d_Curve& AReferenceCurve, |
| 507 | const Standard_Real Tolerance, |
| 508 | const TColStd_Array1OfReal& Parameters, |
| 509 | Standard_Real& MaxDistance) |
| 510 | { |
| 511 | Standard_Integer ii ; |
| 512 | Standard_Real max_squared = 0.0e0, |
| 513 | tolerance_squared = Tolerance * Tolerance, |
| 514 | other_parameter, |
| 515 | para_tolerance, |
| 516 | local_distance_squared ; |
| 517 | gp_Pnt Point1 ; |
| 518 | gp_Pnt Point2 ; |
| 519 | |
| 520 | |
| 521 | |
| 522 | para_tolerance = |
| 523 | AReferenceCurve.Resolution(Tolerance) ; |
| 524 | other_parameter = Parameters(Parameters.Lower()) ; |
| 525 | ACurve.D0(other_parameter, |
| 526 | Point1) ; |
| 527 | Extrema_LocateExtPC a_projector(Point1, |
| 528 | AReferenceCurve, |
| 529 | other_parameter, |
| 530 | para_tolerance) ; |
| 531 | for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) { |
| 532 | ACurve.D0(Parameters(ii), |
| 533 | Point1) ; |
| 534 | AReferenceCurve.D0(Parameters(ii), |
| 535 | Point2) ; |
| 536 | local_distance_squared = |
| 537 | Point1.SquareDistance (Point2) ; |
| 538 | |
| 539 | local_distance_squared = |
| 540 | Point1.SquareDistance (Point2) ; |
| 541 | |
| 542 | |
| 543 | if (local_distance_squared > tolerance_squared) { |
| 544 | |
| 545 | |
| 546 | a_projector.Perform(Point1, |
| 547 | other_parameter) ; |
| 548 | if (a_projector.IsDone()) { |
| 549 | other_parameter = |
| 550 | a_projector.Point().Parameter() ; |
| 551 | AReferenceCurve.D0(other_parameter, |
| 552 | Point2) ; |
| 553 | local_distance_squared = |
| 554 | Point1.SquareDistance (Point2) ; |
| 555 | } |
| 556 | else { |
| 557 | local_distance_squared = 0.0e0 ; |
| 558 | other_parameter = Parameters(ii) ; |
| 559 | } |
| 560 | } |
| 561 | else { |
| 562 | other_parameter = Parameters(ii) ; |
| 563 | } |
| 564 | |
| 565 | |
| 566 | max_squared = Max(max_squared,local_distance_squared) ; |
| 567 | } |
| 568 | if (max_squared > tolerance_squared) { |
| 569 | MaxDistance = sqrt(max_squared) ; |
| 570 | } |
| 571 | else { |
| 572 | MaxDistance = Tolerance ; |
| 573 | } |
| 574 | } |
| 575 | |
| 576 | |
| 577 | |
| 578 | // Aliases: |
| 579 | |
| 580 | // Global data definitions: |
| 581 | |
| 582 | // Methods : |
| 583 | |
| 584 | |
| 585 | //======================================================================= |
| 586 | //function : To3d |
| 587 | //purpose : |
| 588 | //======================================================================= |
| 589 | |
| 590 | Handle(Geom_Curve) GeomLib::To3d (const gp_Ax2& Position, |
| 591 | const Handle(Geom2d_Curve)& Curve2d ) { |
| 592 | Handle(Geom_Curve) Curve3d; |
| 593 | Handle(Standard_Type) KindOfCurve = Curve2d->DynamicType(); |
| 594 | |
| 595 | if (KindOfCurve == STANDARD_TYPE (Geom2d_TrimmedCurve)) { |
| 596 | Handle(Geom2d_TrimmedCurve) Ct = |
| 597 | Handle(Geom2d_TrimmedCurve)::DownCast(Curve2d); |
| 598 | Standard_Real U1 = Ct->FirstParameter (); |
| 599 | Standard_Real U2 = Ct->LastParameter (); |
| 600 | Handle(Geom2d_Curve) CBasis2d = Ct->BasisCurve(); |
| 601 | Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d); |
| 602 | Curve3d = new Geom_TrimmedCurve (CC, U1, U2); |
| 603 | } |
| 604 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_OffsetCurve)) { |
| 605 | Handle(Geom2d_OffsetCurve) Co = |
| 606 | Handle(Geom2d_OffsetCurve)::DownCast(Curve2d); |
| 607 | Standard_Real Offset = Co->Offset(); |
| 608 | Handle(Geom2d_Curve) CBasis2d = Co->BasisCurve(); |
| 609 | Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d); |
| 610 | Curve3d = new Geom_OffsetCurve (CC, Offset, Position.Direction()); |
| 611 | } |
| 612 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_BezierCurve)) { |
| 613 | Handle(Geom2d_BezierCurve) CBez2d = |
| 614 | Handle(Geom2d_BezierCurve)::DownCast (Curve2d); |
| 615 | Standard_Integer Nbpoles = CBez2d->NbPoles (); |
| 616 | TColgp_Array1OfPnt2d Poles2d (1, Nbpoles); |
| 617 | CBez2d->Poles (Poles2d); |
| 618 | TColgp_Array1OfPnt Poles3d (1, Nbpoles); |
| 619 | for (Standard_Integer i = 1; i <= Nbpoles; i++) { |
| 620 | Poles3d (i) = ElCLib::To3d (Position, Poles2d (i)); |
| 621 | } |
| 622 | Handle(Geom_BezierCurve) CBez3d; |
| 623 | if (CBez2d->IsRational()) { |
| 624 | TColStd_Array1OfReal TheWeights (1, Nbpoles); |
| 625 | CBez2d->Weights (TheWeights); |
| 626 | CBez3d = new Geom_BezierCurve (Poles3d, TheWeights); |
| 627 | } |
| 628 | else { |
| 629 | CBez3d = new Geom_BezierCurve (Poles3d); |
| 630 | } |
| 631 | Curve3d = CBez3d; |
| 632 | } |
| 633 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_BSplineCurve)) { |
| 634 | Handle(Geom2d_BSplineCurve) CBSpl2d = |
| 635 | Handle(Geom2d_BSplineCurve)::DownCast (Curve2d); |
| 636 | Standard_Integer Nbpoles = CBSpl2d->NbPoles (); |
| 637 | Standard_Integer Nbknots = CBSpl2d->NbKnots (); |
| 638 | Standard_Integer TheDegree = CBSpl2d->Degree (); |
| 639 | Standard_Boolean IsPeriodic = CBSpl2d->IsPeriodic(); |
| 640 | TColgp_Array1OfPnt2d Poles2d (1, Nbpoles); |
| 641 | CBSpl2d->Poles (Poles2d); |
| 642 | TColgp_Array1OfPnt Poles3d (1, Nbpoles); |
| 643 | for (Standard_Integer i = 1; i <= Nbpoles; i++) { |
| 644 | Poles3d (i) = ElCLib::To3d (Position, Poles2d (i)); |
| 645 | } |
| 646 | TColStd_Array1OfReal TheKnots (1, Nbknots); |
| 647 | TColStd_Array1OfInteger TheMults (1, Nbknots); |
| 648 | CBSpl2d->Knots (TheKnots); |
| 649 | CBSpl2d->Multiplicities (TheMults); |
| 650 | Handle(Geom_BSplineCurve) CBSpl3d; |
| 651 | if (CBSpl2d->IsRational()) { |
| 652 | TColStd_Array1OfReal TheWeights (1, Nbpoles); |
| 653 | CBSpl2d->Weights (TheWeights); |
| 654 | CBSpl3d = new Geom_BSplineCurve (Poles3d, TheWeights, TheKnots, TheMults, TheDegree, IsPeriodic); |
| 655 | } |
| 656 | else { |
| 657 | CBSpl3d = new Geom_BSplineCurve (Poles3d, TheKnots, TheMults, TheDegree, IsPeriodic); |
| 658 | } |
| 659 | Curve3d = CBSpl3d; |
| 660 | } |
| 661 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_Line)) { |
| 662 | Handle(Geom2d_Line) Line2d = Handle(Geom2d_Line)::DownCast (Curve2d); |
| 663 | gp_Lin2d L2d = Line2d->Lin2d(); |
| 664 | gp_Lin L3d = ElCLib::To3d (Position, L2d); |
| 665 | Handle(Geom_Line) GeomL3d = new Geom_Line (L3d); |
| 666 | Curve3d = GeomL3d; |
| 667 | } |
| 668 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_Circle)) { |
| 669 | Handle(Geom2d_Circle) Circle2d = |
| 670 | Handle(Geom2d_Circle)::DownCast (Curve2d); |
| 671 | gp_Circ2d C2d = Circle2d->Circ2d(); |
| 672 | gp_Circ C3d = ElCLib::To3d (Position, C2d); |
| 673 | Handle(Geom_Circle) GeomC3d = new Geom_Circle (C3d); |
| 674 | Curve3d = GeomC3d; |
| 675 | } |
| 676 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_Ellipse)) { |
| 677 | Handle(Geom2d_Ellipse) Ellipse2d = |
| 678 | Handle(Geom2d_Ellipse)::DownCast (Curve2d); |
| 679 | gp_Elips2d E2d = Ellipse2d->Elips2d (); |
| 680 | gp_Elips E3d = ElCLib::To3d (Position, E2d); |
| 681 | Handle(Geom_Ellipse) GeomE3d = new Geom_Ellipse (E3d); |
| 682 | Curve3d = GeomE3d; |
| 683 | } |
| 684 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_Parabola)) { |
| 685 | Handle(Geom2d_Parabola) Parabola2d = |
| 686 | Handle(Geom2d_Parabola)::DownCast (Curve2d); |
| 687 | gp_Parab2d Prb2d = Parabola2d->Parab2d (); |
| 688 | gp_Parab Prb3d = ElCLib::To3d (Position, Prb2d); |
| 689 | Handle(Geom_Parabola) GeomPrb3d = new Geom_Parabola (Prb3d); |
| 690 | Curve3d = GeomPrb3d; |
| 691 | } |
| 692 | else if (KindOfCurve == STANDARD_TYPE (Geom2d_Hyperbola)) { |
| 693 | Handle(Geom2d_Hyperbola) Hyperbola2d = |
| 694 | Handle(Geom2d_Hyperbola)::DownCast (Curve2d); |
| 695 | gp_Hypr2d H2d = Hyperbola2d->Hypr2d (); |
| 696 | gp_Hypr H3d = ElCLib::To3d (Position, H2d); |
| 697 | Handle(Geom_Hyperbola) GeomH3d = new Geom_Hyperbola (H3d); |
| 698 | Curve3d = GeomH3d; |
| 699 | } |
| 700 | else { |
| 701 | Standard_NotImplemented::Raise(); |
| 702 | } |
| 703 | |
| 704 | return Curve3d; |
| 705 | } |
| 706 | |
| 707 | |
| 708 | |
| 709 | //======================================================================= |
| 710 | //function : GTransform |
| 711 | //purpose : |
| 712 | //======================================================================= |
| 713 | |
| 714 | Handle(Geom2d_Curve) GeomLib::GTransform(const Handle(Geom2d_Curve)& Curve, |
| 715 | const gp_GTrsf2d& GTrsf) |
| 716 | { |
| 717 | gp_TrsfForm Form = GTrsf.Form(); |
| 718 | |
| 719 | if ( Form != gp_Other) { |
| 720 | |
| 721 | // Alors, la GTrsf est en fait une Trsf. |
| 722 | // La geometrie des courbes sera alors inchangee. |
| 723 | |
| 724 | Handle(Geom2d_Curve) C = |
| 725 | Handle(Geom2d_Curve)::DownCast(Curve->Transformed(GTrsf.Trsf2d())); |
| 726 | return C; |
| 727 | } |
| 728 | else { |
| 729 | |
| 730 | // Alors, la GTrsf est une other Transformation. |
| 731 | // La geometrie des courbes est alors changee, et les conics devront |
| 732 | // etre converties en BSplines. |
| 733 | |
| 734 | Handle(Standard_Type) TheType = Curve->DynamicType(); |
| 735 | |
| 736 | if ( TheType == STANDARD_TYPE(Geom2d_TrimmedCurve)) { |
| 737 | |
| 738 | // On va recurer sur la BasisCurve |
| 739 | |
| 740 | Handle(Geom2d_TrimmedCurve) C = |
| 741 | Handle(Geom2d_TrimmedCurve)::DownCast(Curve->Copy()); |
| 742 | |
| 743 | Handle(Standard_Type) TheBasisType = (C->BasisCurve())->DynamicType(); |
| 744 | |
| 745 | if (TheBasisType == STANDARD_TYPE(Geom2d_BSplineCurve) || |
| 746 | TheBasisType == STANDARD_TYPE(Geom2d_BezierCurve) ) { |
| 747 | |
| 748 | // Dans ces cas le parametrage est conserve sur la courbe transformee |
| 749 | // on peut donc la trimmer avec les parametres de la courbe de base. |
| 750 | |
| 751 | Standard_Real U1 = C->FirstParameter(); |
| 752 | Standard_Real U2 = C->LastParameter(); |
| 753 | |
| 754 | Handle(Geom2d_TrimmedCurve) result = |
| 755 | new Geom2d_TrimmedCurve(GTransform(C->BasisCurve(), GTrsf), U1,U2); |
| 756 | return result; |
| 757 | } |
| 758 | else if ( TheBasisType == STANDARD_TYPE(Geom2d_Line)) { |
| 759 | |
| 760 | // Dans ce cas, le parametrage n`est plus conserve. |
| 761 | // Il faut recalculer les parametres de Trimming sur la courbe |
| 762 | // resultante. ( Calcul par projection ( ElCLib) des points debut |
| 763 | // et fin transformes) |
| 764 | |
| 765 | Handle(Geom2d_Line) L = |
| 766 | Handle(Geom2d_Line)::DownCast(GTransform(C->BasisCurve(), GTrsf)); |
| 767 | gp_Lin2d Lin = L->Lin2d(); |
| 768 | |
| 769 | gp_Pnt2d P1 = C->StartPoint(); |
| 770 | gp_Pnt2d P2 = C->EndPoint(); |
| 771 | P1.SetXY(GTrsf.Transformed(P1.XY())); |
| 772 | P2.SetXY(GTrsf.Transformed(P2.XY())); |
| 773 | Standard_Real U1 = ElCLib::Parameter(Lin,P1); |
| 774 | Standard_Real U2 = ElCLib::Parameter(Lin,P2); |
| 775 | |
| 776 | Handle(Geom2d_TrimmedCurve) result = |
| 777 | new Geom2d_TrimmedCurve(L,U1,U2); |
| 778 | return result; |
| 779 | } |
| 780 | else if (TheBasisType == STANDARD_TYPE(Geom2d_Circle) || |
| 781 | TheBasisType == STANDARD_TYPE(Geom2d_Ellipse) || |
| 782 | TheBasisType == STANDARD_TYPE(Geom2d_Parabola) || |
| 783 | TheBasisType == STANDARD_TYPE(Geom2d_Hyperbola) ) { |
| 784 | |
| 785 | // Dans ces cas, la geometrie de la courbe n`est pas conservee |
| 786 | // on la convertir en BSpline avant de lui appliquer la Trsf. |
| 787 | |
| 788 | Handle(Geom2d_BSplineCurve) BS = |
| 789 | Geom2dConvert::CurveToBSplineCurve(C); |
| 790 | return GTransform(BS,GTrsf); |
| 791 | } |
| 792 | else { |
| 793 | |
| 794 | // La transformee d`une OffsetCurve vaut ????? Sais pas faire !! |
| 795 | |
| 796 | Handle(Geom2d_Curve) dummy; |
| 797 | return dummy; |
| 798 | } |
| 799 | } |
| 800 | else if ( TheType == STANDARD_TYPE(Geom2d_Line)) { |
| 801 | |
| 802 | Handle(Geom2d_Line) L = |
| 803 | Handle(Geom2d_Line)::DownCast(Curve->Copy()); |
| 804 | gp_Lin2d Lin = L->Lin2d(); |
| 805 | gp_Pnt2d P = Lin.Location(); |
| 806 | gp_Pnt2d PP = L->Value(10.); // pourquoi pas !! |
| 807 | P.SetXY(GTrsf.Transformed(P.XY())); |
| 808 | PP.SetXY(GTrsf.Transformed(PP.XY())); |
| 809 | L->SetLocation(P); |
| 810 | gp_Vec2d V(P,PP); |
| 811 | L->SetDirection(gp_Dir2d(V)); |
| 812 | return L; |
| 813 | } |
| 814 | else if ( TheType == STANDARD_TYPE(Geom2d_BezierCurve)) { |
| 815 | |
| 816 | // Les GTrsf etant des operation lineaires, la transformee d`une courbe |
| 817 | // a poles est la courbe dont les poles sont la transformee des poles |
| 818 | // de la courbe de base. |
| 819 | |
| 820 | Handle(Geom2d_BezierCurve) C = |
| 821 | Handle(Geom2d_BezierCurve)::DownCast(Curve->Copy()); |
| 822 | Standard_Integer NbPoles = C->NbPoles(); |
| 823 | TColgp_Array1OfPnt2d Poles(1,NbPoles); |
| 824 | C->Poles(Poles); |
| 825 | for ( Standard_Integer i = 1; i <= NbPoles; i++) { |
| 826 | Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY())); |
| 827 | C->SetPole(i,Poles(i)); |
| 828 | } |
| 829 | return C; |
| 830 | } |
| 831 | else if ( TheType == STANDARD_TYPE(Geom2d_BSplineCurve)) { |
| 832 | |
| 833 | // Voir commentaire pour les Bezier. |
| 834 | |
| 835 | Handle(Geom2d_BSplineCurve) C = |
| 836 | Handle(Geom2d_BSplineCurve)::DownCast(Curve->Copy()); |
| 837 | Standard_Integer NbPoles = C->NbPoles(); |
| 838 | TColgp_Array1OfPnt2d Poles(1,NbPoles); |
| 839 | C->Poles(Poles); |
| 840 | for ( Standard_Integer i = 1; i <= NbPoles; i++) { |
| 841 | Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY())); |
| 842 | C->SetPole(i,Poles(i)); |
| 843 | } |
| 844 | return C; |
| 845 | } |
| 846 | else if ( TheType == STANDARD_TYPE(Geom2d_Circle) || |
| 847 | TheType == STANDARD_TYPE(Geom2d_Ellipse) ) { |
| 848 | |
| 849 | // Dans ces cas, la geometrie de la courbe n`est pas conservee |
| 850 | // on la convertir en BSpline avant de lui appliquer la Trsf. |
| 851 | |
| 852 | Handle(Geom2d_BSplineCurve) C = |
| 853 | Geom2dConvert::CurveToBSplineCurve(Curve); |
| 854 | return GTransform(C, GTrsf); |
| 855 | } |
| 856 | else if ( TheType == STANDARD_TYPE(Geom2d_Parabola) || |
| 857 | TheType == STANDARD_TYPE(Geom2d_Hyperbola) || |
| 858 | TheType == STANDARD_TYPE(Geom2d_OffsetCurve) ) { |
| 859 | |
| 860 | // On ne sait pas faire : return a null Handle; |
| 861 | |
| 862 | Handle(Geom2d_Curve) dummy; |
| 863 | return dummy; |
| 864 | } |
| 865 | } |
| 866 | |
| 867 | Handle(Geom2d_Curve) WNT__; // portage Windows. |
| 868 | return WNT__; |
| 869 | } |
| 870 | |
| 871 | |
| 872 | //======================================================================= |
| 873 | //function : SameRange |
| 874 | //purpose : |
| 875 | //======================================================================= |
| 876 | void GeomLib::SameRange(const Standard_Real Tolerance, |
| 877 | const Handle(Geom2d_Curve)& CurvePtr, |
| 878 | const Standard_Real FirstOnCurve, |
| 879 | const Standard_Real LastOnCurve, |
| 880 | const Standard_Real RequestedFirst, |
| 881 | const Standard_Real RequestedLast, |
| 882 | Handle(Geom2d_Curve)& NewCurvePtr) |
| 883 | { |
| 884 | if(CurvePtr.IsNull()) Standard_Failure::Raise(); |
| 885 | if (Abs(LastOnCurve - RequestedLast) <= Tolerance && |
| 886 | Abs(FirstOnCurve - RequestedFirst) <= Tolerance) { |
| 887 | NewCurvePtr = CurvePtr; |
| 888 | return; |
| 889 | } |
| 890 | |
| 891 | // the parametrisation lentgh must at least be the same. |
| 892 | if (Abs(LastOnCurve - FirstOnCurve - RequestedLast + RequestedFirst) |
| 893 | <= Tolerance) { |
| 894 | if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Line))) { |
| 895 | Handle(Geom2d_Line) Line = |
| 896 | Handle(Geom2d_Line)::DownCast(CurvePtr->Copy()); |
| 897 | Standard_Real dU = FirstOnCurve - RequestedFirst; |
| 898 | gp_Dir2d D = Line->Direction() ; |
| 899 | Line->Translate(dU * gp_Vec2d(D)); |
| 900 | NewCurvePtr = Line; |
| 901 | } |
| 902 | else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Circle))) { |
| 903 | gp_Trsf2d Trsf; |
| 904 | NewCurvePtr = Handle(Geom2d_Curve)::DownCast(CurvePtr->Copy()); |
| 905 | Handle(Geom2d_Circle) Circ = |
| 906 | Handle(Geom2d_Circle)::DownCast(NewCurvePtr); |
| 907 | gp_Pnt2d P = Circ->Location(); |
| 908 | Standard_Real dU; |
| 909 | if (Circ->Circ2d().IsDirect()) { |
| 910 | dU = FirstOnCurve - RequestedFirst; |
| 911 | } |
| 912 | else { |
| 913 | dU = RequestedFirst - FirstOnCurve; |
| 914 | } |
| 915 | Trsf.SetRotation(P,dU); |
| 916 | NewCurvePtr->Transform(Trsf) ; |
| 917 | } |
| 918 | else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) { |
| 919 | Handle(Geom2d_TrimmedCurve) TC = |
| 920 | Handle(Geom2d_TrimmedCurve)::DownCast(CurvePtr); |
| 921 | GeomLib::SameRange(Tolerance, |
| 922 | TC->BasisCurve(), |
| 923 | FirstOnCurve , LastOnCurve, |
| 924 | RequestedFirst, RequestedLast, |
| 925 | NewCurvePtr); |
| 926 | NewCurvePtr = new Geom2d_TrimmedCurve( NewCurvePtr, RequestedFirst, RequestedLast ); |
| 927 | } |
| 928 | // |
| 929 | // attention a des problemes de limitation : utiliser le MEME test que dans |
| 930 | // Geom2d_TrimmedCurve::SetTrim car sinon comme on risque de relimite sur |
| 931 | // RequestedFirst et RequestedLast on aura un probleme |
| 932 | // |
| 933 | // |
| 934 | else if (Abs(LastOnCurve - FirstOnCurve) > Precision::PConfusion() || |
| 935 | Abs(RequestedLast + RequestedFirst) > Precision::PConfusion()) { |
| 936 | |
| 937 | Handle(Geom2d_TrimmedCurve) TC = |
| 938 | new Geom2d_TrimmedCurve(CurvePtr,FirstOnCurve,LastOnCurve); |
| 939 | |
| 940 | Handle(Geom2d_BSplineCurve) BS = |
| 941 | Geom2dConvert::CurveToBSplineCurve(TC); |
| 942 | TColStd_Array1OfReal Knots(1,BS->NbKnots()); |
| 943 | BS->Knots(Knots); |
| 944 | |
| 945 | BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots); |
| 946 | |
| 947 | BS->SetKnots(Knots); |
| 948 | NewCurvePtr = BS; |
| 949 | } |
| 950 | |
| 951 | } |
| 952 | else { // On segmente le resultat |
| 953 | Handle(Geom2d_TrimmedCurve) TC = |
| 954 | new Geom2d_TrimmedCurve( CurvePtr, FirstOnCurve, LastOnCurve ); |
| 955 | // |
| 956 | Handle(Geom2d_BSplineCurve) BS = |
| 957 | Geom2dConvert::CurveToBSplineCurve(TC); |
| 958 | TColStd_Array1OfReal Knots(1,BS->NbKnots()); |
| 959 | BS->Knots(Knots); |
| 960 | |
| 961 | BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots); |
| 962 | |
| 963 | BS->SetKnots(Knots); |
| 964 | NewCurvePtr = BS; |
| 965 | } |
| 966 | } |
| 967 | |
| 968 | //======================================================================= |
| 969 | //class : GeomLib_CurveOnSurfaceEvaluator |
| 970 | //purpose: The evaluator for the Curve 3D building |
| 971 | //======================================================================= |
| 972 | |
| 973 | class GeomLib_CurveOnSurfaceEvaluator : public AdvApprox_EvaluatorFunction |
| 974 | { |
| 975 | public: |
| 976 | GeomLib_CurveOnSurfaceEvaluator (Adaptor3d_CurveOnSurface& theCurveOnSurface, |
| 977 | Standard_Real theFirst, Standard_Real theLast) |
| 978 | : CurveOnSurface(theCurveOnSurface), FirstParam(theFirst), LastParam(theLast) {} |
| 979 | |
| 980 | virtual void Evaluate (Standard_Integer *Dimension, |
| 981 | Standard_Real StartEnd[2], |
| 982 | Standard_Real *Parameter, |
| 983 | Standard_Integer *DerivativeRequest, |
| 984 | Standard_Real *Result, // [Dimension] |
| 985 | Standard_Integer *ErrorCode); |
| 986 | |
| 987 | private: |
| 988 | Adaptor3d_CurveOnSurface& CurveOnSurface; |
| 989 | Standard_Real FirstParam; |
| 990 | Standard_Real LastParam; |
| 991 | |
| 992 | Handle(Adaptor3d_HCurve) TrimCurve; |
| 993 | }; |
| 994 | |
| 995 | void GeomLib_CurveOnSurfaceEvaluator::Evaluate (Standard_Integer *,/*Dimension*/ |
| 996 | Standard_Real DebutFin[2], |
| 997 | Standard_Real *Parameter, |
| 998 | Standard_Integer *DerivativeRequest, |
| 999 | Standard_Real *Result,// [Dimension] |
| 1000 | Standard_Integer *ReturnCode) |
| 1001 | { |
| 1002 | register Standard_Integer ii ; |
| 1003 | gp_Pnt Point ; |
| 1004 | |
| 1005 | //Gestion des positionnements gauche / droite |
| 1006 | if ((DebutFin[0] != FirstParam) || (DebutFin[1] != LastParam)) |
| 1007 | { |
| 1008 | TrimCurve = CurveOnSurface.Trim(DebutFin[0], DebutFin[1], Precision::PConfusion()); |
| 1009 | FirstParam = DebutFin[0]; |
| 1010 | LastParam = DebutFin[1]; |
| 1011 | } |
| 1012 | |
| 1013 | //Positionemment |
| 1014 | if (*DerivativeRequest == 0) |
| 1015 | { |
| 1016 | TrimCurve->D0((*Parameter), Point) ; |
| 1017 | |
| 1018 | for (ii = 0 ; ii < 3 ; ii++) |
| 1019 | Result[ii] = Point.Coord(ii + 1); |
| 1020 | } |
| 1021 | if (*DerivativeRequest == 1) |
| 1022 | { |
| 1023 | gp_Vec Vector; |
| 1024 | TrimCurve->D1((*Parameter), Point, Vector); |
| 1025 | for (ii = 0 ; ii < 3 ; ii++) |
| 1026 | Result[ii] = Vector.Coord(ii + 1) ; |
| 1027 | } |
| 1028 | if (*DerivativeRequest == 2) |
| 1029 | { |
| 1030 | gp_Vec Vector, VecBis; |
| 1031 | TrimCurve->D2((*Parameter), Point, VecBis, Vector); |
| 1032 | for (ii = 0 ; ii < 3 ; ii++) |
| 1033 | Result[ii] = Vector.Coord(ii + 1) ; |
| 1034 | } |
| 1035 | ReturnCode[0] = 0; |
| 1036 | } |
| 1037 | |
| 1038 | //======================================================================= |
| 1039 | //function : BuildCurve3d |
| 1040 | //purpose : |
| 1041 | //======================================================================= |
| 1042 | |
| 1043 | void GeomLib::BuildCurve3d(const Standard_Real Tolerance, |
| 1044 | Adaptor3d_CurveOnSurface& Curve, |
| 1045 | const Standard_Real FirstParameter, |
| 1046 | const Standard_Real LastParameter, |
| 1047 | Handle(Geom_Curve)& NewCurvePtr, |
| 1048 | Standard_Real& MaxDeviation, |
| 1049 | Standard_Real& AverageDeviation, |
| 1050 | const GeomAbs_Shape Continuity, |
| 1051 | const Standard_Integer MaxDegree, |
| 1052 | const Standard_Integer MaxSegment) |
| 1053 | |
| 1054 | { |
| 1055 | |
| 1056 | |
| 1057 | Standard_Integer curve_not_computed = 1 ; |
| 1058 | MaxDeviation = 0.0e0 ; |
| 1059 | AverageDeviation = 0.0e0 ; |
| 1060 | const Handle(GeomAdaptor_HSurface) & geom_adaptor_surface_ptr = |
| 1061 | Handle(GeomAdaptor_HSurface)::DownCast(Curve.GetSurface()) ; |
| 1062 | const Handle(Geom2dAdaptor_HCurve) & geom_adaptor_curve_ptr = |
| 1063 | Handle(Geom2dAdaptor_HCurve)::DownCast(Curve.GetCurve()) ; |
| 1064 | |
| 1065 | if (! geom_adaptor_curve_ptr.IsNull() && |
| 1066 | ! geom_adaptor_surface_ptr.IsNull()) { |
| 1067 | Handle(Geom_Plane) P ; |
| 1068 | const GeomAdaptor_Surface & geom_surface = |
| 1069 | * (GeomAdaptor_Surface *) &geom_adaptor_surface_ptr->Surface() ; |
| 1070 | |
| 1071 | Handle(Geom_RectangularTrimmedSurface) RT = |
| 1072 | Handle(Geom_RectangularTrimmedSurface):: |
| 1073 | DownCast(geom_surface.Surface()); |
| 1074 | if ( RT.IsNull()) { |
| 1075 | P = Handle(Geom_Plane)::DownCast(geom_surface.Surface()); |
| 1076 | } |
| 1077 | else { |
| 1078 | P = Handle(Geom_Plane)::DownCast(RT->BasisSurface()); |
| 1079 | } |
| 1080 | |
| 1081 | |
| 1082 | if (! P.IsNull()) { |
| 1083 | // compute the 3d curve |
| 1084 | gp_Ax2 axes = P->Position().Ax2(); |
| 1085 | const Geom2dAdaptor_Curve & geom2d_curve = |
| 1086 | * (Geom2dAdaptor_Curve *) & geom_adaptor_curve_ptr->Curve2d() ; |
| 1087 | NewCurvePtr = |
| 1088 | GeomLib::To3d(axes, |
| 1089 | geom2d_curve.Curve()); |
| 1090 | curve_not_computed = 0 ; |
| 1091 | |
| 1092 | } |
| 1093 | } |
| 1094 | if (curve_not_computed) { |
| 1095 | |
| 1096 | // |
| 1097 | // Entree |
| 1098 | // |
| 1099 | Handle(TColStd_HArray1OfReal) Tolerance1DPtr,Tolerance2DPtr; |
| 1100 | Handle(TColStd_HArray1OfReal) Tolerance3DPtr = |
| 1101 | new TColStd_HArray1OfReal(1,1) ; |
| 1102 | Tolerance3DPtr->SetValue(1,Tolerance); |
| 1103 | |
| 1104 | // Recherche des discontinuitees |
| 1105 | Standard_Integer NbIntervalC2 = Curve.NbIntervals(GeomAbs_C2); |
| 1106 | TColStd_Array1OfReal Param_de_decoupeC2 (1, NbIntervalC2+1); |
| 1107 | Curve.Intervals(Param_de_decoupeC2, GeomAbs_C2); |
| 1108 | |
| 1109 | Standard_Integer NbIntervalC3 = Curve.NbIntervals(GeomAbs_C3); |
| 1110 | TColStd_Array1OfReal Param_de_decoupeC3 (1, NbIntervalC3+1); |
| 1111 | Curve.Intervals(Param_de_decoupeC3, GeomAbs_C3); |
| 1112 | |
| 1113 | // Note extension of the parameteric range |
| 1114 | // Pour forcer le Trim au premier appel de l'evaluateur |
| 1115 | GeomLib_CurveOnSurfaceEvaluator ev (Curve, FirstParameter - 1., LastParameter + 1.); |
| 1116 | |
| 1117 | // Approximation avec decoupe preferentiel |
| 1118 | AdvApprox_PrefAndRec Preferentiel(Param_de_decoupeC2, |
| 1119 | Param_de_decoupeC3); |
| 1120 | AdvApprox_ApproxAFunction anApproximator(0, |
| 1121 | 0, |
| 1122 | 1, |
| 1123 | Tolerance1DPtr, |
| 1124 | Tolerance2DPtr, |
| 1125 | Tolerance3DPtr, |
| 1126 | FirstParameter, |
| 1127 | LastParameter, |
| 1128 | Continuity, |
| 1129 | MaxDegree, |
| 1130 | MaxSegment, |
| 1131 | ev, |
| 1132 | // CurveOnSurfaceEvaluator, |
| 1133 | Preferentiel) ; |
| 1134 | |
| 1135 | if (anApproximator.HasResult()) { |
| 1136 | GeomLib_MakeCurvefromApprox |
| 1137 | aCurveBuilder(anApproximator) ; |
| 1138 | |
| 1139 | Handle(Geom_BSplineCurve) aCurvePtr = |
| 1140 | aCurveBuilder.Curve(1) ; |
| 1141 | // On rend les resultats de l'approx |
| 1142 | MaxDeviation = anApproximator.MaxError(3,1) ; |
| 1143 | AverageDeviation = anApproximator.AverageError(3,1) ; |
| 1144 | NewCurvePtr = aCurvePtr ; |
| 1145 | } |
| 1146 | } |
| 1147 | } |
| 1148 | |
| 1149 | //======================================================================= |
| 1150 | //function : AdjustExtremity |
| 1151 | //purpose : |
| 1152 | //======================================================================= |
| 1153 | |
| 1154 | void GeomLib::AdjustExtremity(Handle(Geom_BoundedCurve)& Curve, |
| 1155 | const gp_Pnt& P1, |
| 1156 | const gp_Pnt& P2, |
| 1157 | const gp_Vec& T1, |
| 1158 | const gp_Vec& T2) |
| 1159 | { |
| 1160 | // il faut Convertir l'entree (en preservant si possible le parametrage) |
| 1161 | Handle(Geom_BSplineCurve) aIn, aDef; |
| 1162 | aIn = GeomConvert::CurveToBSplineCurve(Curve, Convert_QuasiAngular); |
| 1163 | |
| 1164 | Standard_Integer ii, jj; |
| 1165 | gp_Pnt P; |
| 1166 | gp_Vec V, Vtan, DV; |
| 1167 | TColgp_Array1OfPnt PolesDef(1,4), Coeffs(1,4); |
| 1168 | TColStd_Array1OfReal FK(1, 8); |
| 1169 | TColStd_Array1OfReal Ti(1, 4); |
| 1170 | TColStd_Array1OfInteger Contact(1, 4); |
| 1171 | |
| 1172 | Ti(1) = Ti(2) = aIn->FirstParameter(); |
| 1173 | Ti(3) = Ti(4) = aIn->LastParameter(); |
| 1174 | Contact(1) = Contact(3) = 0; |
| 1175 | Contact(2) = Contact(4) = 1; |
| 1176 | for (ii=1; ii<=4; ii++) { |
| 1177 | FK(ii) = aIn->FirstParameter(); |
| 1178 | FK(ii) = aIn->LastParameter(); |
| 1179 | } |
| 1180 | |
| 1181 | // Calculs des contraintes de deformations |
| 1182 | aIn->D1(Ti(1), P, V); |
| 1183 | PolesDef(1).ChangeCoord() = P1.XYZ()-P.XYZ(); |
| 1184 | Vtan = T1; |
| 1185 | Vtan.Normalize(); |
| 1186 | DV = Vtan * (Vtan * V) - V; |
| 1187 | PolesDef(2).ChangeCoord() = (Ti(4)-Ti(1))*DV.XYZ(); |
| 1188 | |
| 1189 | aIn->D1(Ti(4), P, V); |
| 1190 | PolesDef(3).ChangeCoord() = P2.XYZ()-P.XYZ(); |
| 1191 | Vtan = T2; |
| 1192 | Vtan.Normalize(); |
| 1193 | DV = Vtan * (Vtan * V) - V; |
| 1194 | PolesDef(4).ChangeCoord() = (Ti(4)-Ti(1))* DV.XYZ(); |
| 1195 | |
| 1196 | // Interpolation des contraintes |
| 1197 | math_Matrix Mat(1, 4, 1, 4); |
| 1198 | if (!PLib::HermiteCoefficients(0., 1., 1, 1, Mat)) |
| 1199 | Standard_ConstructionError::Raise(); |
| 1200 | |
| 1201 | for (jj=1; jj<=4; jj++) { |
| 1202 | gp_XYZ aux(0.,0.,0.); |
| 1203 | for (ii=1; ii<=4; ii++) { |
| 1204 | aux.SetLinearForm(Mat(ii,jj), PolesDef(ii).XYZ(), aux); |
| 1205 | } |
| 1206 | Coeffs(jj).SetXYZ(aux); |
| 1207 | } |
| 1208 | |
| 1209 | PLib::CoefficientsPoles(Coeffs, PLib::NoWeights(), |
| 1210 | PolesDef, PLib::NoWeights()); |
| 1211 | |
| 1212 | // Ajout de la deformation |
| 1213 | TColStd_Array1OfReal K(1, 2); |
| 1214 | TColStd_Array1OfInteger M(1, 2); |
| 1215 | K(1) = Ti(1); |
| 1216 | K(2) = Ti(4); |
| 1217 | M.Init(4); |
| 1218 | |
| 1219 | aDef = new (Geom_BSplineCurve) (PolesDef, K, M, 3); |
| 1220 | if (aIn->Degree() < 3) aIn->IncreaseDegree(3); |
| 1221 | else aDef->IncreaseDegree(aIn->Degree()); |
| 1222 | |
| 1223 | for (ii=2; ii<aIn->NbKnots(); ii++) { |
| 1224 | aDef->InsertKnot(aIn->Knot(ii), aIn->Multiplicity(ii)); |
| 1225 | } |
| 1226 | |
| 1227 | if (aDef->NbPoles() != aIn->NbPoles()) |
| 1228 | Standard_ConstructionError::Raise("Inconsistent poles's number"); |
| 1229 | |
| 1230 | for (ii=1; ii<=aDef->NbPoles(); ii++) { |
| 1231 | P = aIn->Pole(ii); |
| 1232 | P.ChangeCoord() += aDef->Pole(ii).XYZ(); |
| 1233 | aIn->SetPole(ii, P); |
| 1234 | } |
| 1235 | Curve = aIn; |
| 1236 | } |
| 1237 | //======================================================================= |
| 1238 | //function : ExtendCurveToPoint |
| 1239 | //purpose : |
| 1240 | //======================================================================= |
| 1241 | |
| 1242 | void GeomLib::ExtendCurveToPoint(Handle(Geom_BoundedCurve)& Curve, |
| 1243 | const gp_Pnt& Point, |
| 1244 | const Standard_Integer Continuity, |
| 1245 | const Standard_Boolean After) |
| 1246 | { |
| 1247 | if(Continuity < 1 || Continuity > 3) return; |
| 1248 | Standard_Integer size = Continuity + 2; |
| 1249 | Standard_Real Ubord, Tol=1.e-6; |
| 1250 | math_Matrix MatCoefs(1,size, 1,size); |
| 1251 | Standard_Real Lambda, L1; |
| 1252 | Standard_Integer ii, jj; |
| 1253 | gp_Vec d1, d2, d3; |
| 1254 | gp_Pnt p0; |
| 1255 | // il faut Convertir l'entree (en preservant si possible le parametrage) |
| 1256 | GeomConvert_CompCurveToBSplineCurve Concat(Curve, Convert_QuasiAngular); |
| 1257 | |
| 1258 | // Les contraintes de constructions |
| 1259 | TColgp_Array1OfXYZ Cont(1,size); |
| 1260 | if (After) { |
| 1261 | Ubord = Curve->LastParameter(); |
| 1262 | |
| 1263 | } |
| 1264 | else { |
| 1265 | Ubord = Curve->FirstParameter(); |
| 1266 | } |
| 1267 | PLib::HermiteCoefficients(0, 1, // Les Bornes |
| 1268 | Continuity, 0, // Les Ordres de contraintes |
| 1269 | MatCoefs); |
| 1270 | |
| 1271 | Curve->D3(Ubord, p0, d1, d2, d3); |
| 1272 | if (!After) { // Inversion du parametrage |
| 1273 | d1 *= -1; |
| 1274 | d3 *= -1; |
| 1275 | } |
| 1276 | |
| 1277 | L1 = p0.Distance(Point); |
| 1278 | if (L1 > Tol) { |
| 1279 | // Lambda est le ratio qu'il faut appliquer a la derive de la courbe |
| 1280 | // pour obtenir la derive du prolongement (fixe arbitrairement a la |
| 1281 | // longueur du segment bout de la courbe - point cible. |
| 1282 | // On essai d'avoir sur le prolongement la vitesse moyenne que l'on |
| 1283 | // a sur la courbe. |
| 1284 | gp_Vec daux; |
| 1285 | gp_Pnt pp; |
| 1286 | Standard_Real f= Curve->FirstParameter(), t, dt, norm; |
| 1287 | dt = (Curve->LastParameter()-f)/9; |
| 1288 | norm = d1.Magnitude(); |
| 1289 | for (ii=1, t=f+dt; ii<=8; ii++, t+=dt) { |
| 1290 | Curve->D1(t, pp, daux); |
| 1291 | norm += daux.Magnitude(); |
| 1292 | } |
| 1293 | norm /= 9; |
| 1294 | dt = d1.Magnitude() / norm; |
| 1295 | if ((dt<1.5) && (dt>0.75)) { // Le bord est dans la moyenne on le garde |
| 1296 | Lambda = ((Standard_Real)1) / Max (d1.Magnitude() / L1, Tol); |
| 1297 | } |
| 1298 | else { |
| 1299 | Lambda = ((Standard_Real)1) / Max (norm / L1, Tol); |
| 1300 | } |
| 1301 | } |
| 1302 | else { |
| 1303 | return; // Pas d'extension |
| 1304 | } |
| 1305 | |
| 1306 | // Optimisation du Lambda |
| 1307 | math_Matrix Cons(1, 3, 1, size); |
| 1308 | Cons(1,1) = p0.X(); Cons(2,1) = p0.Y(); Cons(3,1) = p0.Z(); |
| 1309 | Cons(1,2) = d1.X(); Cons(2,2) = d1.Y(); Cons(3,2) = d1.Z(); |
| 1310 | Cons(1,size) = Point.X(); Cons(2,size) = Point.Y(); Cons(3,size) = Point.Z(); |
| 1311 | if (Continuity >= 2) { |
| 1312 | Cons(1,3) = d2.X(); Cons(2,3) = d2.Y(); Cons(3,3) = d2.Z(); |
| 1313 | } |
| 1314 | if (Continuity >= 3) { |
| 1315 | Cons(1,4) = d3.X(); Cons(2,4) = d3.Y(); Cons(3,4) = d3.Z(); |
| 1316 | } |
| 1317 | ComputeLambda(Cons, MatCoefs, L1, Lambda); |
| 1318 | |
| 1319 | // Construction dans la Base Polynomiale |
| 1320 | Cont(1) = p0.XYZ(); |
| 1321 | Cont(2) = d1.XYZ() * Lambda; |
| 1322 | if(Continuity >= 2) Cont(3) = d2.XYZ() * Pow(Lambda,2); |
| 1323 | if(Continuity >= 3) Cont(4) = d3.XYZ() * Pow(Lambda,3); |
| 1324 | Cont(size) = Point.XYZ(); |
| 1325 | |
| 1326 | |
| 1327 | TColgp_Array1OfPnt ExtrapPoles(1, size); |
| 1328 | TColgp_Array1OfPnt ExtraCoeffs(1, size); |
| 1329 | |
| 1330 | gp_Pnt PNull(0.,0.,0.); |
| 1331 | ExtraCoeffs.Init(PNull); |
| 1332 | for (ii=1; ii<=size; ii++) { |
| 1333 | for (jj=1; jj<=size; jj++) { |
| 1334 | ExtraCoeffs(jj).ChangeCoord() += MatCoefs(ii,jj)*Cont(ii); |
| 1335 | } |
| 1336 | } |
| 1337 | |
| 1338 | // Convertion Dans la Base de Bernstein |
| 1339 | PLib::CoefficientsPoles(ExtraCoeffs, PLib::NoWeights(), |
| 1340 | ExtrapPoles, PLib::NoWeights()); |
| 1341 | |
| 1342 | Handle(Geom_BezierCurve) Bezier = new (Geom_BezierCurve) (ExtrapPoles); |
| 1343 | |
| 1344 | Standard_Real dist = ExtrapPoles(1).Distance(p0); |
| 1345 | Standard_Boolean Ok; |
| 1346 | Tol += dist; |
| 1347 | |
| 1348 | // Concatenation |
| 1349 | Ok = Concat.Add(Bezier, Tol, After); |
| 1350 | if (!Ok) Standard_ConstructionError::Raise("ExtendCurveToPoint"); |
| 1351 | |
| 1352 | Curve = Concat.BSplineCurve(); |
| 1353 | } |
| 1354 | |
| 1355 | |
| 1356 | //======================================================================= |
| 1357 | //function : ExtendKPart |
| 1358 | //purpose : Extension par longueur des surfaces cannonique |
| 1359 | //======================================================================= |
| 1360 | static Standard_Boolean |
| 1361 | ExtendKPart(Handle(Geom_RectangularTrimmedSurface)& Surface, |
| 1362 | const Standard_Real Length, |
| 1363 | const Standard_Boolean InU, |
| 1364 | const Standard_Boolean After) |
| 1365 | { |
| 1366 | |
| 1367 | if (Surface.IsNull()) return Standard_False; |
| 1368 | |
| 1369 | Standard_Boolean Ok=Standard_True; |
| 1370 | Standard_Real Uf, Ul, Vf, Vl; |
| 1371 | Handle(Geom_Surface) Support = Surface->BasisSurface(); |
| 1372 | GeomAbs_SurfaceType Type; |
| 1373 | |
| 1374 | Surface->Bounds(Uf, Ul, Vf, Vl); |
| 1375 | GeomAdaptor_Surface AS(Surface); |
| 1376 | Type = AS.GetType(); |
| 1377 | |
| 1378 | if (InU) { |
| 1379 | switch(Type) { |
| 1380 | case GeomAbs_Plane : |
| 1381 | { |
| 1382 | if (After) Ul+=Length; |
| 1383 | else Uf-=Length; |
| 1384 | Surface = new (Geom_RectangularTrimmedSurface) |
| 1385 | (Support, Uf, Ul, Vf, Vl); |
| 1386 | break; |
| 1387 | } |
| 1388 | |
| 1389 | default: |
| 1390 | Ok = Standard_False; |
| 1391 | } |
| 1392 | } |
| 1393 | else { |
| 1394 | switch(Type) { |
| 1395 | case GeomAbs_Plane : |
| 1396 | case GeomAbs_Cylinder : |
| 1397 | case GeomAbs_SurfaceOfExtrusion : |
| 1398 | { |
| 1399 | if (After) Vl+=Length; |
| 1400 | else Vf-=Length; |
| 1401 | Surface = new (Geom_RectangularTrimmedSurface) |
| 1402 | (Support, Uf, Ul, Vf, Vl); |
| 1403 | break; |
| 1404 | } |
| 1405 | default: |
| 1406 | Ok = Standard_False; |
| 1407 | } |
| 1408 | } |
| 1409 | |
| 1410 | return Ok; |
| 1411 | } |
| 1412 | |
| 1413 | //======================================================================= |
| 1414 | //function : ExtendSurfByLength |
| 1415 | //purpose : |
| 1416 | //======================================================================= |
| 1417 | void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface, |
| 1418 | const Standard_Real Length, |
| 1419 | const Standard_Integer Continuity, |
| 1420 | const Standard_Boolean InU, |
| 1421 | const Standard_Boolean After) |
| 1422 | { |
| 1423 | if(Continuity < 0 || Continuity > 3) return; |
| 1424 | Standard_Integer Cont = Continuity; |
| 1425 | |
| 1426 | // Kpart ? |
| 1427 | Handle(Geom_RectangularTrimmedSurface) TS = |
| 1428 | Handle(Geom_RectangularTrimmedSurface)::DownCast (Surface); |
| 1429 | if (ExtendKPart(TS,Length, InU, After) ) { |
| 1430 | Surface = TS; |
| 1431 | return; |
| 1432 | } |
| 1433 | |
| 1434 | // format BSplineSurface avec un degre suffisant pour la continuite voulue |
| 1435 | Handle(Geom_BSplineSurface) BS = |
| 1436 | Handle(Geom_BSplineSurface)::DownCast (Surface); |
| 1437 | if (BS.IsNull()) { |
| 1438 | //BS = GeomConvert::SurfaceToBSplineSurface(Surface); |
| 1439 | Standard_Real Tol = Precision::Confusion(); //1.e-4; |
| 1440 | GeomAbs_Shape UCont = GeomAbs_C1, VCont = GeomAbs_C1; |
| 1441 | Standard_Integer degU = 14, degV = 14; |
| 1442 | Standard_Integer nmax = 16; |
| 1443 | Standard_Integer thePrec = 1; |
| 1444 | GeomConvert_ApproxSurface theApprox(Surface,Tol,UCont,VCont,degU,degV,nmax,thePrec); |
| 1445 | if (theApprox.HasResult()) |
| 1446 | BS = theApprox.Surface(); |
| 1447 | else |
| 1448 | BS = GeomConvert::SurfaceToBSplineSurface(Surface); |
| 1449 | } |
| 1450 | if (InU&&(BS->UDegree()<Continuity+1)) |
| 1451 | BS->IncreaseDegree(Continuity+1,BS->VDegree()); |
| 1452 | if (!InU&&(BS->VDegree()<Continuity+1)) |
| 1453 | BS->IncreaseDegree(BS->UDegree(),Continuity+1); |
| 1454 | |
| 1455 | // si BS etait periodique dans le sens de l'extension, elle ne le sera plus |
| 1456 | if ( (InU&&(BS->IsUPeriodic())) || (!InU&&(BS->IsVPeriodic())) ) { |
| 1457 | Standard_Real U0,U1,V0,V1; |
| 1458 | BS->Bounds(U0,U1,V0,V1); |
| 1459 | BS->Segment(U0,U1,V0,V1); |
| 1460 | } |
| 1461 | |
| 1462 | |
| 1463 | // IFV Fix OCC bug 0022694 - wrong result extrapolating rational surfaces |
| 1464 | // Standard_Boolean rational = ( InU && BS->IsURational() ) |
| 1465 | // || ( !InU && BS->IsVRational() ) ; |
| 1466 | Standard_Boolean rational = (BS->IsURational() || BS->IsVRational()); |
| 1467 | Standard_Boolean NullWeight; |
| 1468 | Standard_Real EpsW = 10*Precision::PConfusion(); |
| 1469 | Standard_Integer gap = 3; |
| 1470 | if ( rational ) gap++; |
| 1471 | |
| 1472 | |
| 1473 | |
| 1474 | Standard_Integer Cdeg = 0, Cdim = 0, NbP = 0, Ksize = 0, Psize = 1; |
| 1475 | Standard_Integer ii, jj, ipole, Kount; |
| 1476 | Standard_Real Tbord, lambmin=Length; |
| 1477 | Standard_Real * Padr = NULL; |
| 1478 | Standard_Boolean Ok; |
| 1479 | Handle(TColStd_HArray1OfReal) FKnots, Point, lambda, Tgte, Poles; |
| 1480 | |
| 1481 | |
| 1482 | |
| 1483 | |
| 1484 | for (Kount=0, Ok=Standard_False; Kount<=2 && !Ok; Kount++) { |
| 1485 | // transformation de la surface en une BSpline non rationnelle a une variable |
| 1486 | // de degre UDegree ou VDegree et de dimension 3 ou 4 x NbVpoles ou NbUpoles |
| 1487 | // le nombre de poles egal a NbUpoles ou NbVpoles |
| 1488 | // ATTENTION : dans le cas rationnel, un point de coordonnees (x,y,z) |
| 1489 | // et de poids w devient un point de coordonnees (wx, wy, wz, w ) |
| 1490 | |
| 1491 | |
| 1492 | if (InU) { |
| 1493 | Cdeg = BS->UDegree(); |
| 1494 | NbP = BS->NbUPoles(); |
| 1495 | Cdim = BS->NbVPoles() * gap; |
| 1496 | } |
| 1497 | else { |
| 1498 | Cdeg = BS->VDegree(); |
| 1499 | NbP = BS->NbVPoles(); |
| 1500 | Cdim = BS->NbUPoles() * gap; |
| 1501 | } |
| 1502 | |
| 1503 | // les noeuds plats |
| 1504 | Ksize = NbP + Cdeg + 1; |
| 1505 | FKnots = new (TColStd_HArray1OfReal) (1,Ksize); |
| 1506 | if (InU) |
| 1507 | BS->UKnotSequence(FKnots->ChangeArray1()); |
| 1508 | else |
| 1509 | BS->VKnotSequence(FKnots->ChangeArray1()); |
| 1510 | |
| 1511 | // le parametre du noeud de raccord |
| 1512 | if (After) |
| 1513 | Tbord = FKnots->Value(FKnots->Upper()-Cdeg); |
| 1514 | else |
| 1515 | Tbord = FKnots->Value(FKnots->Lower()+Cdeg); |
| 1516 | |
| 1517 | // les poles |
| 1518 | Psize = Cdim * NbP; |
| 1519 | Poles = new (TColStd_HArray1OfReal) (1,Psize); |
| 1520 | |
| 1521 | if (InU) { |
| 1522 | for (ii=1,ipole=1; ii<=NbP; ii++) { |
| 1523 | for (jj=1;jj<=BS->NbVPoles();jj++) { |
| 1524 | Poles->SetValue(ipole, BS->Pole(ii,jj).X()); |
| 1525 | Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y()); |
| 1526 | Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z()); |
| 1527 | if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj)); |
| 1528 | ipole+=gap; |
| 1529 | } |
| 1530 | } |
| 1531 | } |
| 1532 | else { |
| 1533 | for (jj=1,ipole=1; jj<=NbP; jj++) { |
| 1534 | for (ii=1;ii<=BS->NbUPoles();ii++) { |
| 1535 | Poles->SetValue(ipole, BS->Pole(ii,jj).X()); |
| 1536 | Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y()); |
| 1537 | Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z()); |
| 1538 | if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj)); |
| 1539 | ipole+=gap; |
| 1540 | } |
| 1541 | } |
| 1542 | } |
| 1543 | Padr = (Standard_Real *) &Poles->ChangeValue(1); |
| 1544 | |
| 1545 | // calcul du point de raccord et de la tangente |
| 1546 | Point = new (TColStd_HArray1OfReal)(1,Cdim); |
| 1547 | Tgte = new (TColStd_HArray1OfReal)(1,Cdim); |
| 1548 | lambda = new (TColStd_HArray1OfReal)(1,Cdim); |
| 1549 | |
| 1550 | Standard_Boolean periodic_flag = Standard_False ; |
| 1551 | Standard_Integer extrap_mode[2], derivative_request = Max(Continuity,1); |
| 1552 | extrap_mode[0] = extrap_mode[1] = Cdeg; |
| 1553 | TColStd_Array1OfReal Result(1, Cdim * (derivative_request+1)) ; |
| 1554 | |
| 1555 | TColStd_Array1OfReal& tgte = Tgte->ChangeArray1(); |
| 1556 | TColStd_Array1OfReal& point = Point->ChangeArray1(); |
| 1557 | TColStd_Array1OfReal& lamb = lambda->ChangeArray1(); |
| 1558 | |
| 1559 | Standard_Real * Radr = (Standard_Real *) &Result(1) ; |
| 1560 | |
| 1561 | BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0], |
| 1562 | Cdeg,FKnots->Array1(),Cdim,*Padr,*Radr); |
| 1563 | Ok = Standard_True; |
| 1564 | for (ii=1;ii<=Cdim;ii++) { |
| 1565 | point(ii) = Result(ii); |
| 1566 | tgte(ii) = Result(ii+Cdim); |
| 1567 | } |
| 1568 | |
| 1569 | // calcul de la contrainte a atteindre |
| 1570 | |
| 1571 | gp_Vec CurT, OldT; |
| 1572 | |
| 1573 | Standard_Real NTgte, val, Tgtol = 1.e-12, OldN = 0.0; |
| 1574 | if (rational) { |
| 1575 | for (ii=gap;ii<=Cdim;ii+=gap) { |
| 1576 | tgte(ii) = 0.; |
| 1577 | } |
| 1578 | for (ii=gap;ii<=Cdim;ii+=gap) { |
| 1579 | CurT.SetCoord(tgte(ii-3),tgte(ii-2), tgte(ii-1)); |
| 1580 | NTgte=CurT.Magnitude(); |
| 1581 | if (NTgte>Tgtol) { |
| 1582 | val = Length/NTgte; |
| 1583 | // Attentions aux Cas ou le segment donne par les poles |
| 1584 | // est oppose au sens de la derive |
| 1585 | // Exemple: Certaine portions de tore. |
| 1586 | if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) { |
| 1587 | Ok = Standard_False; |
| 1588 | } |
| 1589 | |
| 1590 | lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = val; |
| 1591 | lamb(ii) = 0.; |
| 1592 | lambmin = Min(lambmin, val); |
| 1593 | } |
| 1594 | else { |
| 1595 | lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = 0.; |
| 1596 | lamb(ii) = 0.; |
| 1597 | } |
| 1598 | OldT = CurT; |
| 1599 | OldN = NTgte; |
| 1600 | } |
| 1601 | } |
| 1602 | else { |
| 1603 | for (ii=gap;ii<=Cdim;ii+=gap) { |
| 1604 | CurT.SetCoord(tgte(ii-2),tgte(ii-1), tgte(ii)); |
| 1605 | NTgte=CurT.Magnitude(); |
| 1606 | if (NTgte>Tgtol) { |
| 1607 | val = Length/NTgte; |
| 1608 | // Attentions aux Cas ou le segment donne par les poles |
| 1609 | // est oppose au sens de la derive |
| 1610 | // Exemple: Certaine portion de tore. |
| 1611 | if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) { |
| 1612 | Ok = Standard_False; |
| 1613 | } |
| 1614 | lamb(ii) = lamb(ii-1) = lamb(ii-2) = val; |
| 1615 | lambmin = Min(lambmin, val); |
| 1616 | } |
| 1617 | else { |
| 1618 | lamb(ii) =lamb(ii-1) = lamb(ii-2) = 0.; |
| 1619 | } |
| 1620 | OldT = CurT; |
| 1621 | OldN = NTgte; |
| 1622 | } |
| 1623 | } |
| 1624 | if (!Ok && Kount<2) { |
| 1625 | // On augmente le degre de l'iso bord afin de rapprocher les poles de la surface |
| 1626 | // Et on ressaye |
| 1627 | if (InU) BS->IncreaseDegree(BS->UDegree(), BS->VDegree()+2); |
| 1628 | else BS->IncreaseDegree(BS->UDegree()+2, BS->VDegree()); |
| 1629 | } |
| 1630 | } |
| 1631 | |
| 1632 | |
| 1633 | TColStd_Array1OfReal ConstraintPoint(1,Cdim); |
| 1634 | if (After) { |
| 1635 | for (ii=1;ii<=Cdim;ii++) { |
| 1636 | ConstraintPoint(ii) = Point->Value(ii) + lambda->Value(ii)*Tgte->Value(ii); |
| 1637 | } |
| 1638 | } |
| 1639 | else { |
| 1640 | for (ii=1;ii<=Cdim;ii++) { |
| 1641 | ConstraintPoint(ii) = Point->Value(ii) - lambda->Value(ii)*Tgte->Value(ii); |
| 1642 | } |
| 1643 | } |
| 1644 | |
| 1645 | // cas particulier du rationnel |
| 1646 | if (rational) { |
| 1647 | for (ipole=1;ipole<=Psize;ipole+=gap) { |
| 1648 | Poles->ChangeValue(ipole) *= Poles->Value(ipole+3); |
| 1649 | Poles->ChangeValue(ipole+1) *= Poles->Value(ipole+3); |
| 1650 | Poles->ChangeValue(ipole+2) *= Poles->Value(ipole+3); |
| 1651 | } |
| 1652 | for (ii=1;ii<=Cdim;ii+=gap) { |
| 1653 | ConstraintPoint(ii) *= ConstraintPoint(ii+3); |
| 1654 | ConstraintPoint(ii+1) *= ConstraintPoint(ii+3); |
| 1655 | ConstraintPoint(ii+2) *= ConstraintPoint(ii+3); |
| 1656 | } |
| 1657 | } |
| 1658 | |
| 1659 | // tableaux necessaires pour l'extension |
| 1660 | Standard_Integer Ksize2 = Ksize+Cdeg, NbPoles, NbKnots = 0; |
| 1661 | TColStd_Array1OfReal FK(1, Ksize2) ; |
| 1662 | Standard_Real * FKRadr = &FK(1); |
| 1663 | |
| 1664 | Standard_Integer Psize2 = Psize+Cdeg*Cdim; |
| 1665 | TColStd_Array1OfReal PRes(1, Psize2) ; |
| 1666 | Standard_Real * PRadr = &PRes(1); |
| 1667 | Standard_Real ww; |
| 1668 | Standard_Boolean ExtOk = Standard_False; |
| 1669 | Handle(TColgp_HArray2OfPnt) NewPoles; |
| 1670 | Handle(TColStd_HArray2OfReal) NewWeights; |
| 1671 | |
| 1672 | |
| 1673 | for (Kount=1; Kount<=5 && !ExtOk; Kount++) { |
| 1674 | // extension |
| 1675 | BSplCLib::TangExtendToConstraint(FKnots->Array1(), |
| 1676 | lambmin,NbP,*Padr, |
| 1677 | Cdim,Cdeg, |
| 1678 | ConstraintPoint, Cont, After, |
| 1679 | NbPoles, NbKnots,*FKRadr, *PRadr); |
| 1680 | |
| 1681 | // recopie des poles du resultat sous forme de points 3D et de poids |
| 1682 | Standard_Integer NU, NV, indice ; |
| 1683 | if (InU) { |
| 1684 | NU = NbPoles; |
| 1685 | NV = BS->NbVPoles(); |
| 1686 | } |
| 1687 | else { |
| 1688 | NU = BS->NbUPoles(); |
| 1689 | NV = NbPoles; |
| 1690 | } |
| 1691 | |
| 1692 | NewPoles = new (TColgp_HArray2OfPnt)(1,NU,1,NV); |
| 1693 | TColgp_Array2OfPnt& NewP = NewPoles->ChangeArray2(); |
| 1694 | NewWeights = new (TColStd_HArray2OfReal) (1,NU,1,NV); |
| 1695 | TColStd_Array2OfReal& NewW = NewWeights->ChangeArray2(); |
| 1696 | |
| 1697 | if (!rational) NewW.Init(1.); |
| 1698 | NullWeight= Standard_False; |
| 1699 | |
| 1700 | if (InU) { |
| 1701 | for (ii=1; ii<=NU && !NullWeight; ii++) { |
| 1702 | for (jj=1; jj<=NV && !NullWeight; jj++) { |
| 1703 | indice = 1+(ii-1)*Cdim+(jj-1)*gap; |
| 1704 | NewP(ii,jj).SetCoord(1,PRes(indice)); |
| 1705 | NewP(ii,jj).SetCoord(2,PRes(indice+1)); |
| 1706 | NewP(ii,jj).SetCoord(3,PRes(indice+2)); |
| 1707 | if (rational) { |
| 1708 | ww = PRes(indice+3); |
| 1709 | if (ww < EpsW) { |
| 1710 | NullWeight = Standard_True; |
| 1711 | } |
| 1712 | else { |
| 1713 | NewW(ii,jj) = ww; |
| 1714 | NewP(ii,jj).ChangeCoord() /= ww; |
| 1715 | } |
| 1716 | } |
| 1717 | } |
| 1718 | } |
| 1719 | } |
| 1720 | else { |
| 1721 | for (jj=1; jj<=NV && !NullWeight; jj++) { |
| 1722 | for (ii=1; ii<=NU && !NullWeight; ii++) { |
| 1723 | indice = 1+(ii-1)*gap+(jj-1)*Cdim; |
| 1724 | NewP(ii,jj).SetCoord(1,PRes(indice)); |
| 1725 | NewP(ii,jj).SetCoord(2,PRes(indice+1)); |
| 1726 | NewP(ii,jj).SetCoord(3,PRes(indice+2)); |
| 1727 | if (rational) { |
| 1728 | ww = PRes(indice+3); |
| 1729 | if (ww < EpsW) { |
| 1730 | NullWeight = Standard_True; |
| 1731 | } |
| 1732 | else { |
| 1733 | NewW(ii,jj) = ww; |
| 1734 | NewP(ii,jj).ChangeCoord() /= ww; |
| 1735 | } |
| 1736 | } |
| 1737 | } |
| 1738 | } |
| 1739 | } |
| 1740 | |
| 1741 | if (NullWeight) { |
| 1742 | #ifdef OCCT_DEBUG |
| 1743 | cout << "Echec de l'Extension rationnelle" << endl; |
| 1744 | #endif |
| 1745 | lambmin /= 3.; |
| 1746 | NullWeight = Standard_False; |
| 1747 | } |
| 1748 | else { |
| 1749 | ExtOk = Standard_True; |
| 1750 | } |
| 1751 | } |
| 1752 | |
| 1753 | |
| 1754 | // recopie des noeuds plats sous forme de noeuds avec leurs multiplicites |
| 1755 | // calcul des degres du resultat |
| 1756 | Standard_Integer Usize = BS->NbUKnots(), Vsize = BS->NbVKnots(), UDeg, VDeg; |
| 1757 | if (InU) |
| 1758 | Usize++; |
| 1759 | else |
| 1760 | Vsize++; |
| 1761 | TColStd_Array1OfReal UKnots(1,Usize); |
| 1762 | TColStd_Array1OfReal VKnots(1,Vsize); |
| 1763 | TColStd_Array1OfInteger UMults(1,Usize); |
| 1764 | TColStd_Array1OfInteger VMults(1,Vsize); |
| 1765 | TColStd_Array1OfReal FKRes(1, NbKnots); |
| 1766 | |
| 1767 | for (ii=1; ii<=NbKnots; ii++) |
| 1768 | FKRes(ii) = FK(ii); |
| 1769 | |
| 1770 | if (InU) { |
| 1771 | BSplCLib::Knots(FKRes, UKnots, UMults); |
| 1772 | UDeg = Cdeg; |
| 1773 | UMults(Usize) = UDeg+1; // Petite verrue utile quand la continuite |
| 1774 | // n'est pas ok. |
| 1775 | BS->VKnots(VKnots); |
| 1776 | BS->VMultiplicities(VMults); |
| 1777 | VDeg = BS->VDegree(); |
| 1778 | } |
| 1779 | else { |
| 1780 | BSplCLib::Knots(FKRes, VKnots, VMults); |
| 1781 | VDeg = Cdeg; |
| 1782 | VMults(Vsize) = VDeg+1; |
| 1783 | BS->UKnots(UKnots); |
| 1784 | BS->UMultiplicities(UMults); |
| 1785 | UDeg = BS->UDegree(); |
| 1786 | } |
| 1787 | |
| 1788 | // construction de la surface BSpline resultat |
| 1789 | Handle(Geom_BSplineSurface) Res = |
| 1790 | new (Geom_BSplineSurface) (NewPoles->Array2(), |
| 1791 | NewWeights->Array2(), |
| 1792 | UKnots,VKnots, |
| 1793 | UMults,VMults, |
| 1794 | UDeg,VDeg, |
| 1795 | BS->IsUPeriodic(), |
| 1796 | BS->IsVPeriodic()); |
| 1797 | Surface = Res; |
| 1798 | } |
| 1799 | |
| 1800 | //======================================================================= |
| 1801 | //function : Inertia |
| 1802 | //purpose : |
| 1803 | //======================================================================= |
| 1804 | void GeomLib::Inertia(const TColgp_Array1OfPnt& Points, |
| 1805 | gp_Pnt& Bary, |
| 1806 | gp_Dir& XDir, |
| 1807 | gp_Dir& YDir, |
| 1808 | Standard_Real& Xgap, |
| 1809 | Standard_Real& Ygap, |
| 1810 | Standard_Real& Zgap) |
| 1811 | { |
| 1812 | gp_XYZ GB(0., 0., 0.), Diff; |
| 1813 | // gp_Vec A,B,C,D; |
| 1814 | |
| 1815 | Standard_Integer i,nb=Points.Length(); |
| 1816 | GB.SetCoord(0.,0.,0.); |
| 1817 | for (i=1; i<=nb; i++) |
| 1818 | GB += Points(i).XYZ(); |
| 1819 | |
| 1820 | GB /= nb; |
| 1821 | |
| 1822 | math_Matrix M (1, 3, 1, 3); |
| 1823 | M.Init(0.); |
| 1824 | for (i=1; i<=nb; i++) { |
| 1825 | Diff.SetLinearForm(-1, Points(i).XYZ(), GB); |
| 1826 | M(1,1) += Diff.X() * Diff.X(); |
| 1827 | M(2,2) += Diff.Y() * Diff.Y(); |
| 1828 | M(3,3) += Diff.Z() * Diff.Z(); |
| 1829 | M(1,2) += Diff.X() * Diff.Y(); |
| 1830 | M(1,3) += Diff.X() * Diff.Z(); |
| 1831 | M(2,3) += Diff.Y() * Diff.Z(); |
| 1832 | } |
| 1833 | |
| 1834 | M(2,1)=M(1,2) ; |
| 1835 | M(3,1)=M(1,3) ; |
| 1836 | M(3,2)=M(2,3) ; |
| 1837 | |
| 1838 | M /= nb; |
| 1839 | |
| 1840 | math_Jacobi J(M); |
| 1841 | if (!J.IsDone()) { |
| 1842 | #ifdef OCCT_DEBUG |
| 1843 | cout << "Erreur dans Jacobbi" << endl; |
| 1844 | M.Dump(cout); |
| 1845 | #endif |
| 1846 | } |
| 1847 | |
| 1848 | Standard_Real n1,n2,n3; |
| 1849 | |
| 1850 | n1=J.Value(1); |
| 1851 | n2=J.Value(2); |
| 1852 | n3=J.Value(3); |
| 1853 | |
| 1854 | Standard_Real r1 = Min(Min(n1,n2),n3), r2; |
| 1855 | Standard_Integer m1, m2, m3; |
| 1856 | if (r1==n1) { |
| 1857 | m1 = 1; |
| 1858 | r2 = Min(n2,n3); |
| 1859 | if (r2==n2) { |
| 1860 | m2 = 2; |
| 1861 | m3 = 3; |
| 1862 | } |
| 1863 | else { |
| 1864 | m2 = 3; |
| 1865 | m3 = 2; |
| 1866 | } |
| 1867 | } |
| 1868 | else { |
| 1869 | if (r1==n2) { |
| 1870 | m1 = 2 ; |
| 1871 | r2 = Min(n1,n3); |
| 1872 | if (r2==n1) { |
| 1873 | m2 = 1; |
| 1874 | m3 = 3; |
| 1875 | } |
| 1876 | else { |
| 1877 | m2 = 3; |
| 1878 | m3 = 1; |
| 1879 | } |
| 1880 | } |
| 1881 | else { |
| 1882 | m1 = 3 ; |
| 1883 | r2 = Min(n1,n2); |
| 1884 | if (r2==n1) { |
| 1885 | m2 = 1; |
| 1886 | m3 = 2; |
| 1887 | } |
| 1888 | else { |
| 1889 | m2 = 2; |
| 1890 | m3 = 1; |
| 1891 | } |
| 1892 | } |
| 1893 | } |
| 1894 | |
| 1895 | math_Vector V2(1,3),V3(1,3); |
| 1896 | J.Vector(m2,V2); |
| 1897 | J.Vector(m3,V3); |
| 1898 | |
| 1899 | Bary.SetXYZ(GB); |
| 1900 | XDir.SetCoord(V3(1),V3(2),V3(3)); |
| 1901 | YDir.SetCoord(V2(1),V2(2),V2(3)); |
| 1902 | |
| 1903 | Zgap = sqrt(Abs(J.Value(m1))); |
| 1904 | Ygap = sqrt(Abs(J.Value(m2))); |
| 1905 | Xgap = sqrt(Abs(J.Value(m3))); |
| 1906 | } |
| 1907 | //======================================================================= |
| 1908 | //function : AxeOfInertia |
| 1909 | //purpose : |
| 1910 | //======================================================================= |
| 1911 | void GeomLib::AxeOfInertia(const TColgp_Array1OfPnt& Points, |
| 1912 | gp_Ax2& Axe, |
| 1913 | Standard_Boolean& IsSingular, |
| 1914 | const Standard_Real Tol) |
| 1915 | { |
| 1916 | gp_Pnt Bary; |
| 1917 | gp_Dir OX,OY,OZ; |
| 1918 | Standard_Real gx, gy, gz; |
| 1919 | |
| 1920 | GeomLib::Inertia(Points, Bary, OX, OY, gx, gy, gz); |
| 1921 | |
| 1922 | if (gy*Points.Length()<=Tol) { |
| 1923 | gp_Ax2 axe (Bary, OX); |
| 1924 | OY = axe.XDirection(); |
| 1925 | IsSingular = Standard_True; |
| 1926 | } |
| 1927 | else { |
| 1928 | IsSingular = Standard_False; |
| 1929 | } |
| 1930 | |
| 1931 | OZ = OX^OY; |
| 1932 | gp_Ax2 TheAxe(Bary, OZ, OX); |
| 1933 | Axe = TheAxe; |
| 1934 | } |
| 1935 | |
| 1936 | //======================================================================= |
| 1937 | //function : CanBeTreated |
| 1938 | //purpose : indicates if the surface can be treated(if the conditions are |
| 1939 | // filled) and need to be treated(if the surface hasn't been yet |
| 1940 | // treated or if the surface is rationnal and non periodic) |
| 1941 | //======================================================================= |
| 1942 | |
| 1943 | static Standard_Boolean CanBeTreated(Handle(Geom_BSplineSurface)& BSurf) |
| 1944 | |
| 1945 | {Standard_Integer i; |
| 1946 | Standard_Real lambda; //proportionnality coefficient |
| 1947 | Standard_Boolean AlreadyTreated=Standard_True; |
| 1948 | |
| 1949 | if (!BSurf->IsURational()||(BSurf->IsUPeriodic())) |
| 1950 | return Standard_False; |
| 1951 | else { |
| 1952 | lambda=(BSurf->Weight(1,1)/BSurf->Weight(BSurf->NbUPoles(),1)); |
| 1953 | for (i=1;i<=BSurf->NbVPoles();i++) //test of the proportionnality of the denominator on the boundaries |
| 1954 | if ((BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))<(1-Precision::Confusion()))|| |
| 1955 | (BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))>(1+Precision::Confusion()))) |
| 1956 | return Standard_False; |
| 1957 | i=1; |
| 1958 | while ((AlreadyTreated) && (i<=BSurf->NbVPoles())){ //tests if the surface has already been treated |
| 1959 | if (((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))<(1-Precision::Confusion()))|| |
| 1960 | ((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))>(1+Precision::Confusion()))|| |
| 1961 | ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))<(1-Precision::Confusion()))|| |
| 1962 | ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))>(1+Precision::Confusion()))) |
| 1963 | AlreadyTreated=Standard_False; |
| 1964 | i++; |
| 1965 | } |
| 1966 | if (AlreadyTreated) |
| 1967 | return Standard_False; |
| 1968 | } |
| 1969 | return Standard_True; |
| 1970 | } |
| 1971 | |
| 1972 | //======================================================================= |
| 1973 | //class : law_evaluator |
| 1974 | //purpose : usefull to estimate the value of a function of 2 variables |
| 1975 | //======================================================================= |
| 1976 | |
| 1977 | class law_evaluator : public BSplSLib_EvaluatorFunction |
| 1978 | { |
| 1979 | |
| 1980 | public: |
| 1981 | |
| 1982 | law_evaluator (const GeomLib_DenominatorMultiplierPtr theDenominatorPtr) |
| 1983 | : myDenominator (theDenominatorPtr) {} |
| 1984 | |
| 1985 | virtual void Evaluate (const Standard_Integer theDerivativeRequest, |
| 1986 | const Standard_Real theUParameter, |
| 1987 | const Standard_Real theVParameter, |
| 1988 | Standard_Real& theResult, |
| 1989 | Standard_Integer& theErrorCode) const |
| 1990 | { |
| 1991 | if ((myDenominator != NULL) && (theDerivativeRequest == 0)) |
| 1992 | { |
| 1993 | theResult = myDenominator->Value (theUParameter, theVParameter); |
| 1994 | theErrorCode = 0; |
| 1995 | } |
| 1996 | else |
| 1997 | { |
| 1998 | theErrorCode = 1; |
| 1999 | } |
| 2000 | } |
| 2001 | |
| 2002 | private: |
| 2003 | |
| 2004 | GeomLib_DenominatorMultiplierPtr myDenominator; |
| 2005 | |
| 2006 | }; |
| 2007 | |
| 2008 | //======================================================================= |
| 2009 | //function : CheckIfKnotExists |
| 2010 | //purpose : true if the knot already exists in the knot sequence |
| 2011 | //======================================================================= |
| 2012 | |
| 2013 | static Standard_Boolean CheckIfKnotExists(const TColStd_Array1OfReal& surface_knots, |
| 2014 | const Standard_Real knot) |
| 2015 | |
| 2016 | {Standard_Integer i; |
| 2017 | for (i=1;i<=surface_knots.Length();i++) |
| 2018 | if ((surface_knots(i)-Precision::Confusion()<=knot)&&(surface_knots(i)+Precision::Confusion()>=knot)) |
| 2019 | return Standard_True; |
| 2020 | return Standard_False; |
| 2021 | } |
| 2022 | |
| 2023 | //======================================================================= |
| 2024 | //function : AddAKnot |
| 2025 | //purpose : add a knot and its multiplicity to the knot sequence. This knot |
| 2026 | // will be C2 and the degree is increased of deltasurface_degree |
| 2027 | //======================================================================= |
| 2028 | |
| 2029 | static void AddAKnot(const TColStd_Array1OfReal& knots, |
| 2030 | const TColStd_Array1OfInteger& mults, |
| 2031 | const Standard_Real knotinserted, |
| 2032 | const Standard_Integer deltasurface_degree, |
| 2033 | const Standard_Integer finalsurfacedegree, |
| 2034 | Handle(TColStd_HArray1OfReal) & newknots, |
| 2035 | Handle(TColStd_HArray1OfInteger) & newmults) |
| 2036 | |
| 2037 | {Standard_Integer i; |
| 2038 | |
| 2039 | newknots=new TColStd_HArray1OfReal(1,knots.Length()+1); |
| 2040 | newmults=new TColStd_HArray1OfInteger(1,knots.Length()+1); |
| 2041 | i=1; |
| 2042 | while (knots(i)<knotinserted){ |
| 2043 | newknots->SetValue(i,knots(i)); |
| 2044 | newmults->SetValue(i,mults(i)+deltasurface_degree); |
| 2045 | i++; |
| 2046 | } |
| 2047 | newknots->SetValue(i,knotinserted); //insertion of the new knot |
| 2048 | newmults->SetValue(i,finalsurfacedegree-2); |
| 2049 | i++; |
| 2050 | while (i<=newknots->Length()){ |
| 2051 | newknots->SetValue(i,knots(i-1)); |
| 2052 | newmults->SetValue(i,mults(i-1)+deltasurface_degree); |
| 2053 | i++; |
| 2054 | } |
| 2055 | } |
| 2056 | |
| 2057 | //======================================================================= |
| 2058 | //function : Sort |
| 2059 | //purpose : give the new flat knots(u or v) of the surface |
| 2060 | //======================================================================= |
| 2061 | |
| 2062 | static void BuildFlatKnot(const TColStd_Array1OfReal& surface_knots, |
| 2063 | const TColStd_Array1OfInteger& surface_mults, |
| 2064 | const Standard_Integer deltasurface_degree, |
| 2065 | const Standard_Integer finalsurface_degree, |
| 2066 | const Standard_Real knotmin, |
| 2067 | const Standard_Real knotmax, |
| 2068 | Handle(TColStd_HArray1OfReal)& ResultKnots, |
| 2069 | Handle(TColStd_HArray1OfInteger)& ResultMults) |
| 2070 | |
| 2071 | { |
| 2072 | Standard_Integer i; |
| 2073 | |
| 2074 | if (CheckIfKnotExists(surface_knots,knotmin) && |
| 2075 | CheckIfKnotExists(surface_knots,knotmax)){ |
| 2076 | ResultKnots=new TColStd_HArray1OfReal(1,surface_knots.Length()); |
| 2077 | ResultMults=new TColStd_HArray1OfInteger(1,surface_knots.Length()); |
| 2078 | for (i=1;i<=surface_knots.Length();i++){ |
| 2079 | ResultKnots->SetValue(i,surface_knots(i)); |
| 2080 | ResultMults->SetValue(i,surface_mults(i)+deltasurface_degree); |
| 2081 | } |
| 2082 | } |
| 2083 | else{ |
| 2084 | if ((CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax))) |
| 2085 | AddAKnot(surface_knots,surface_mults,knotmax,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults); |
| 2086 | else{ |
| 2087 | if ((!CheckIfKnotExists(surface_knots,knotmin))&&(CheckIfKnotExists(surface_knots,knotmax))) |
| 2088 | AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults); |
| 2089 | else{ |
| 2090 | if ((!CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax))&& |
| 2091 | (knotmin==knotmax)){ |
| 2092 | AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults); |
| 2093 | } |
| 2094 | else{ |
| 2095 | Handle(TColStd_HArray1OfReal) IntermedKnots; |
| 2096 | Handle(TColStd_HArray1OfInteger) IntermedMults; |
| 2097 | AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,IntermedKnots,IntermedMults); |
| 2098 | AddAKnot(IntermedKnots->ChangeArray1(),IntermedMults->ChangeArray1(),knotmax,0,finalsurface_degree,ResultKnots,ResultMults); |
| 2099 | } |
| 2100 | } |
| 2101 | } |
| 2102 | } |
| 2103 | } |
| 2104 | |
| 2105 | //======================================================================= |
| 2106 | //function : FunctionMultiply |
| 2107 | //purpose : multiply the surface BSurf by a(u,v) (law_evaluator) on its |
| 2108 | // numerator and denominator |
| 2109 | //======================================================================= |
| 2110 | |
| 2111 | static void FunctionMultiply(Handle(Geom_BSplineSurface)& BSurf, |
| 2112 | const Standard_Real knotmin, |
| 2113 | const Standard_Real knotmax) |
| 2114 | |
| 2115 | {TColStd_Array1OfReal surface_u_knots(1,BSurf->NbUKnots()) ; |
| 2116 | TColStd_Array1OfInteger surface_u_mults(1,BSurf->NbUKnots()) ; |
| 2117 | TColStd_Array1OfReal surface_v_knots(1,BSurf->NbVKnots()) ; |
| 2118 | TColStd_Array1OfInteger surface_v_mults(1,BSurf->NbVKnots()) ; |
| 2119 | TColgp_Array2OfPnt surface_poles(1,BSurf->NbUPoles(), |
| 2120 | 1,BSurf->NbVPoles()) ; |
| 2121 | TColStd_Array2OfReal surface_weights(1,BSurf->NbUPoles(), |
| 2122 | 1,BSurf->NbVPoles()) ; |
| 2123 | Standard_Integer i,j,k,status,new_num_u_poles,new_num_v_poles,length=0; |
| 2124 | Handle(TColStd_HArray1OfReal) newuknots,newvknots; |
| 2125 | Handle(TColStd_HArray1OfInteger) newumults,newvmults; |
| 2126 | |
| 2127 | BSurf->UKnots(surface_u_knots) ; |
| 2128 | BSurf->UMultiplicities(surface_u_mults) ; |
| 2129 | BSurf->VKnots(surface_v_knots) ; |
| 2130 | BSurf->VMultiplicities(surface_v_mults) ; |
| 2131 | BSurf->Poles(surface_poles) ; |
| 2132 | BSurf->Weights(surface_weights) ; |
| 2133 | |
| 2134 | TColStd_Array1OfReal Knots(1,2); |
| 2135 | TColStd_Array1OfInteger Mults(1,2); |
| 2136 | Handle(TColStd_HArray1OfReal) NewKnots; |
| 2137 | Handle(TColStd_HArray1OfInteger) NewMults; |
| 2138 | |
| 2139 | Knots(1)=0; |
| 2140 | Knots(2)=1; |
| 2141 | Mults(1)=4; |
| 2142 | Mults(2)=4; |
| 2143 | BuildFlatKnot(Knots,Mults,0,3,knotmin,knotmax,NewKnots,NewMults); |
| 2144 | |
| 2145 | for (i=1;i<=NewMults->Length();i++) |
| 2146 | length+=NewMults->Value(i); |
| 2147 | TColStd_Array1OfReal FlatKnots(1,length); |
| 2148 | BSplCLib::KnotSequence(NewKnots->ChangeArray1(),NewMults->ChangeArray1(),FlatKnots); |
| 2149 | |
| 2150 | GeomLib_DenominatorMultiplier aDenominator (BSurf, FlatKnots); |
| 2151 | |
| 2152 | BuildFlatKnot(surface_u_knots, |
| 2153 | surface_u_mults, |
| 2154 | 3, |
| 2155 | BSurf->UDegree()+3, |
| 2156 | knotmin, |
| 2157 | knotmax, |
| 2158 | newuknots, |
| 2159 | newumults); |
| 2160 | BuildFlatKnot(surface_v_knots, |
| 2161 | surface_v_mults, |
| 2162 | BSurf->VDegree(), |
| 2163 | 2*(BSurf->VDegree()), |
| 2164 | 1.0, |
| 2165 | 0.0, |
| 2166 | newvknots, |
| 2167 | newvmults); |
| 2168 | length=0; |
| 2169 | for (i=1;i<=newumults->Length();i++) |
| 2170 | length+=newumults->Value(i); |
| 2171 | new_num_u_poles=(length-BSurf->UDegree()-3-1); |
| 2172 | TColStd_Array1OfReal newuflatknots(1,length); |
| 2173 | length=0; |
| 2174 | for (i=1;i<=newvmults->Length();i++) |
| 2175 | length+=newvmults->Value(i); |
| 2176 | new_num_v_poles=(length-2*BSurf->VDegree()-1); |
| 2177 | TColStd_Array1OfReal newvflatknots(1,length); |
| 2178 | |
| 2179 | TColgp_Array2OfPnt NewNumerator(1,new_num_u_poles,1,new_num_v_poles); |
| 2180 | TColStd_Array2OfReal NewDenominator(1,new_num_u_poles,1,new_num_v_poles); |
| 2181 | |
| 2182 | BSplCLib::KnotSequence(newuknots->ChangeArray1(),newumults->ChangeArray1(),newuflatknots); |
| 2183 | BSplCLib::KnotSequence(newvknots->ChangeArray1(),newvmults->ChangeArray1(),newvflatknots); |
| 2184 | //POP pour WNT |
| 2185 | law_evaluator ev (&aDenominator); |
| 2186 | // BSplSLib::FunctionMultiply(law_evaluator, //multiplication |
| 2187 | BSplSLib::FunctionMultiply(ev, //multiplication |
| 2188 | BSurf->UDegree(), |
| 2189 | BSurf->VDegree(), |
| 2190 | surface_u_knots, |
| 2191 | surface_v_knots, |
| 2192 | surface_u_mults, |
| 2193 | surface_v_mults, |
| 2194 | surface_poles, |
| 2195 | surface_weights, |
| 2196 | newuflatknots, |
| 2197 | newvflatknots, |
| 2198 | BSurf->UDegree()+3, |
| 2199 | 2*(BSurf->VDegree()), |
| 2200 | NewNumerator, |
| 2201 | NewDenominator, |
| 2202 | status); |
| 2203 | if (status!=0) |
| 2204 | Standard_ConstructionError::Raise("GeomLib Multiplication Error") ; |
| 2205 | for (i = 1 ; i <= new_num_u_poles ; i++) { |
| 2206 | for (j = 1 ; j <= new_num_v_poles ; j++) { |
| 2207 | for (k = 1 ; k <= 3 ; k++) { |
| 2208 | NewNumerator(i,j).SetCoord(k,NewNumerator(i,j).Coord(k)/NewDenominator(i,j)) ; |
| 2209 | } |
| 2210 | } |
| 2211 | } |
| 2212 | BSurf= new Geom_BSplineSurface(NewNumerator, |
| 2213 | NewDenominator, |
| 2214 | newuknots->ChangeArray1(), |
| 2215 | newvknots->ChangeArray1(), |
| 2216 | newumults->ChangeArray1(), |
| 2217 | newvmults->ChangeArray1(), |
| 2218 | BSurf->UDegree()+3, |
| 2219 | 2*(BSurf->VDegree()) ); |
| 2220 | } |
| 2221 | |
| 2222 | //======================================================================= |
| 2223 | //function : CancelDenominatorDerivative1D |
| 2224 | //purpose : cancel the denominator derivative in one direction |
| 2225 | //======================================================================= |
| 2226 | |
| 2227 | static void CancelDenominatorDerivative1D(Handle(Geom_BSplineSurface) & BSurf) |
| 2228 | |
| 2229 | {Standard_Integer i,j; |
| 2230 | Standard_Real uknotmin=1.0,uknotmax=0.0, |
| 2231 | x,y, |
| 2232 | startu_value, |
| 2233 | endu_value; |
| 2234 | TColStd_Array1OfReal BSurf_u_knots(1,BSurf->NbUKnots()) ; |
| 2235 | |
| 2236 | startu_value=BSurf->UKnot(1); |
| 2237 | endu_value=BSurf->UKnot(BSurf->NbUKnots()); |
| 2238 | BSurf->UKnots(BSurf_u_knots) ; |
| 2239 | BSplCLib::Reparametrize(0.0,1.0,BSurf_u_knots); |
| 2240 | BSurf->SetUKnots(BSurf_u_knots); //reparametrisation of the surface |
| 2241 | Handle(Geom_BSplineCurve) BCurve; |
| 2242 | TColStd_Array1OfReal BCurveWeights(1,BSurf->NbUPoles()); |
| 2243 | TColgp_Array1OfPnt BCurvePoles(1,BSurf->NbUPoles()); |
| 2244 | TColStd_Array1OfReal BCurveKnots(1,BSurf->NbUKnots()); |
| 2245 | TColStd_Array1OfInteger BCurveMults(1,BSurf->NbUKnots()); |
| 2246 | |
| 2247 | if (CanBeTreated(BSurf)){ |
| 2248 | for (i=1;i<=BSurf->NbVPoles();i++){ //loop on each pole function |
| 2249 | x=1.0;y=0.0; |
| 2250 | for (j=1;j<=BSurf->NbUPoles();j++){ |
| 2251 | BCurveWeights(j)=BSurf->Weight(j,i); |
| 2252 | BCurvePoles(j)=BSurf->Pole(j,i); |
| 2253 | } |
| 2254 | BSurf->UKnots(BCurveKnots); |
| 2255 | BSurf->UMultiplicities(BCurveMults); |
| 2256 | BCurve = new Geom_BSplineCurve(BCurvePoles, //building of a pole function |
| 2257 | BCurveWeights, |
| 2258 | BCurveKnots, |
| 2259 | BCurveMults, |
| 2260 | BSurf->UDegree()); |
| 2261 | Hermit::Solutionbis(BCurve,x,y,Precision::Confusion(),Precision::Confusion()); |
| 2262 | if (x<uknotmin) |
| 2263 | uknotmin=x; //uknotmin,uknotmax:extremal knots |
| 2264 | if ((x!=1.0)&&(x>uknotmax)) |
| 2265 | uknotmax=x; |
| 2266 | if ((y!=0.0)&&(y<uknotmin)) |
| 2267 | uknotmin=y; |
| 2268 | if (y>uknotmax) |
| 2269 | uknotmax=y; |
| 2270 | } |
| 2271 | |
| 2272 | FunctionMultiply(BSurf,uknotmin,uknotmax); //multiplication |
| 2273 | |
| 2274 | BSurf->UKnots(BSurf_u_knots) ; |
| 2275 | BSplCLib::Reparametrize(startu_value,endu_value,BSurf_u_knots); |
| 2276 | BSurf->SetUKnots(BSurf_u_knots); |
| 2277 | } |
| 2278 | } |
| 2279 | |
| 2280 | //======================================================================= |
| 2281 | //function : CancelDenominatorDerivative |
| 2282 | //purpose : |
| 2283 | //======================================================================= |
| 2284 | |
| 2285 | void GeomLib::CancelDenominatorDerivative(Handle(Geom_BSplineSurface) & BSurf, |
| 2286 | const Standard_Boolean udirection, |
| 2287 | const Standard_Boolean vdirection) |
| 2288 | |
| 2289 | {if (udirection && !vdirection) |
| 2290 | CancelDenominatorDerivative1D(BSurf); |
| 2291 | else{ |
| 2292 | if (!udirection && vdirection) { |
| 2293 | BSurf->ExchangeUV(); |
| 2294 | CancelDenominatorDerivative1D(BSurf); |
| 2295 | BSurf->ExchangeUV(); |
| 2296 | } |
| 2297 | else{ |
| 2298 | if (udirection && vdirection){ //optimize the treatment |
| 2299 | if (BSurf->UDegree()<=BSurf->VDegree()){ |
| 2300 | CancelDenominatorDerivative1D(BSurf); |
| 2301 | BSurf->ExchangeUV(); |
| 2302 | CancelDenominatorDerivative1D(BSurf); |
| 2303 | BSurf->ExchangeUV(); |
| 2304 | } |
| 2305 | else{ |
| 2306 | BSurf->ExchangeUV(); |
| 2307 | CancelDenominatorDerivative1D(BSurf); |
| 2308 | BSurf->ExchangeUV(); |
| 2309 | CancelDenominatorDerivative1D(BSurf); |
| 2310 | } |
| 2311 | } |
| 2312 | } |
| 2313 | } |
| 2314 | } |
| 2315 | |
| 2316 | //======================================================================= |
| 2317 | //function : NormEstim |
| 2318 | //purpose : |
| 2319 | //======================================================================= |
| 2320 | |
| 2321 | Standard_Integer GeomLib::NormEstim(const Handle(Geom_Surface)& S, |
| 2322 | const gp_Pnt2d& UV, |
| 2323 | const Standard_Real Tol, gp_Dir& N) |
| 2324 | { |
| 2325 | gp_Vec DU, DV; |
| 2326 | gp_Pnt DummyPnt; |
| 2327 | Standard_Real aTol2 = Square(Tol); |
| 2328 | |
| 2329 | S->D1(UV.X(), UV.Y(), DummyPnt, DU, DV); |
| 2330 | |
| 2331 | Standard_Real MDU = DU.SquareMagnitude(), MDV = DV.SquareMagnitude(); |
| 2332 | |
| 2333 | if(MDU >= aTol2 && MDV >= aTol2) { |
| 2334 | gp_Vec Norm = DU^DV; |
| 2335 | Standard_Real Magn = Norm.SquareMagnitude(); |
| 2336 | if(Magn < aTol2) return 3; |
| 2337 | |
| 2338 | //Magn = sqrt(Magn); |
| 2339 | N.SetXYZ(Norm.XYZ()); |
| 2340 | |
| 2341 | return 0; |
| 2342 | } |
| 2343 | else { |
| 2344 | gp_Vec D2U, D2V, D2UV; |
| 2345 | Standard_Boolean isDone; |
| 2346 | CSLib_NormalStatus aStatus; |
| 2347 | gp_Dir aNormal; |
| 2348 | |
| 2349 | S->D2(UV.X(), UV.Y(), DummyPnt, DU, DV, D2U, D2V, D2UV); |
| 2350 | CSLib::Normal(DU, DV, D2U, D2V, D2UV, Tol, isDone, aStatus, aNormal); |
| 2351 | |
| 2352 | if (isDone) { |
| 2353 | Standard_Real Umin, Umax, Vmin, Vmax; |
| 2354 | Standard_Real step = 1.0e-5; |
| 2355 | Standard_Real eps = 1.0e-16; |
| 2356 | Standard_Real sign = -1.0; |
| 2357 | |
| 2358 | S->Bounds(Umin, Umax, Vmin, Vmax); |
| 2359 | |
| 2360 | // check for cone apex singularity point |
| 2361 | if ((UV.Y() > Vmin + step) && (UV.Y() < Vmax - step)) |
| 2362 | { |
| 2363 | gp_Dir aNormal1, aNormal2; |
| 2364 | Standard_Real aConeSingularityAngleEps = 1.0e-4; |
| 2365 | S->D1(UV.X(), UV.Y() - sign * step, DummyPnt, DU, DV); |
| 2366 | if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) { |
| 2367 | aNormal1 = DU^DV; |
| 2368 | S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV); |
| 2369 | if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) { |
| 2370 | aNormal2 = DU^DV; |
| 2371 | if (aNormal1.IsOpposite(aNormal2, aConeSingularityAngleEps)) |
| 2372 | return 2; |
| 2373 | } |
| 2374 | } |
| 2375 | } |
| 2376 | |
| 2377 | // Along V |
| 2378 | if(MDU < aTol2 && MDV >= aTol2) { |
| 2379 | if ((Vmax - UV.Y()) > (UV.Y() - Vmin)) |
| 2380 | sign = 1.0; |
| 2381 | S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV); |
| 2382 | gp_Vec Norm = DU^DV; |
| 2383 | if (Norm.SquareMagnitude() < eps) { |
| 2384 | Standard_Real sign1 = -1.0; |
| 2385 | if ((Umax - UV.X()) > (UV.X() - Umin)) |
| 2386 | sign1 = 1.0; |
| 2387 | S->D1(UV.X() + sign1 * step, UV.Y() + sign * step, DummyPnt, DU, DV); |
| 2388 | Norm = DU^DV; |
| 2389 | } |
| 2390 | if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0)) |
| 2391 | aNormal.Reverse(); |
| 2392 | } |
| 2393 | |
| 2394 | // Along U |
| 2395 | if(MDV < aTol2 && MDU >= aTol2) { |
| 2396 | if ((Umax - UV.X()) > (UV.X() - Umin)) |
| 2397 | sign = 1.0; |
| 2398 | S->D1(UV.X() + sign * step, UV.Y(), DummyPnt, DU, DV); |
| 2399 | gp_Vec Norm = DU^DV; |
| 2400 | if (Norm.SquareMagnitude() < eps) { |
| 2401 | Standard_Real sign1 = -1.0; |
| 2402 | if ((Vmax - UV.Y()) > (UV.Y() - Vmin)) |
| 2403 | sign1 = 1.0; |
| 2404 | S->D1(UV.X() + sign * step, UV.Y() + sign1 * step, DummyPnt, DU, DV); |
| 2405 | Norm = DU^DV; |
| 2406 | } |
| 2407 | if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0)) |
| 2408 | aNormal.Reverse(); |
| 2409 | } |
| 2410 | |
| 2411 | // quasysingular |
| 2412 | if ((aStatus == CSLib_D1NuIsNull) || (aStatus == CSLib_D1NvIsNull) || |
| 2413 | (aStatus == CSLib_D1NuIsParallelD1Nv)) { |
| 2414 | N.SetXYZ(aNormal.XYZ()); |
| 2415 | return 1; |
| 2416 | } |
| 2417 | // conical |
| 2418 | if (aStatus == CSLib_InfinityOfSolutions) |
| 2419 | return 2; |
| 2420 | } |
| 2421 | // computation is impossible |
| 2422 | else { |
| 2423 | // conical |
| 2424 | if (aStatus == CSLib_D1NIsNull) { |
| 2425 | return 2; |
| 2426 | } |
| 2427 | return 3; |
| 2428 | } |
| 2429 | } |
| 2430 | return 3; |
| 2431 | } |
| 2432 | |
| 2433 | |