projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023024: Update headers of OCCT files
[occt.git] / src / BRepLib / BRepLib_FindSurface.cxx
CommitLineData
b311480e 1// Created on: 1994-07-22
2// Created by: Remi LEQUETTE
3// Copyright (c) 1994-1999 Matra Datavision
4// Copyright (c) 1999-2012 OPEN CASCADE SAS
5//
6// The content of this file is subject to the Open CASCADE Technology Public
7// License Version 6.5 (the "License"). You may not use the content of this file
8// except in compliance with the License. Please obtain a copy of the License
9// at http://www.opencascade.org and read it completely before using this file.
10//
11// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13//
14// The Original Code and all software distributed under the License is
15// distributed on an "AS IS" basis, without warranty of any kind, and the
16// Initial Developer hereby disclaims all such warranties, including without
17// limitation, any warranties of merchantability, fitness for a particular
18// purpose or non-infringement. Please see the License for the specific terms
19// and conditions governing the rights and limitations under the License.
20
7fd59977 21
22#include <BRepLib_FindSurface.ixx>
23
24#include <Precision.hxx>
25#include <math_Matrix.hxx>
26#include <math_Vector.hxx>
27#include <math_Gauss.hxx>
28
29#include <gp_Lin.hxx>
30#include <gp_Circ.hxx>
31#include <gp_Elips.hxx>
32#include <gp_Hypr.hxx>
33#include <gp_Parab.hxx>
34#include <gp_Ax2.hxx>
35#include <gp_Ax3.hxx>
36#include <gp_Vec.hxx>
37#include <TColgp_SequenceOfPnt.hxx>
38#include <TColStd_SequenceOfReal.hxx>
39#include <TColgp_Array1OfPnt.hxx>
40#include <TColgp_HArray1OfPnt.hxx>
41#include <Geom_Plane.hxx>
42
43#include <TopoDS.hxx>
44#include <TopExp_Explorer.hxx>
45#include <BRep_Tool.hxx>
46#include <BRepAdaptor_Curve.hxx>
47
48#include <GeomLib.hxx>
49#include <Geom2d_Curve.hxx>
50#include <Geom_BezierCurve.hxx>
51#include <Geom_BSplineCurve.hxx>
52
53//=======================================================================
54//function : Controle
55//purpose :
56//=======================================================================
57static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints,
58 const Handle(Geom_Plane)& thePlane)
59{
60 Standard_Real dfMaxDist=0.;
61 Standard_Real a,b,c,d, dist;
62 Standard_Integer ii;
63 thePlane->Coefficients(a,b,c,d);
64 for (ii=1; ii<=thePoints.Length(); ii++) {
65 const gp_XYZ& xyz = thePoints(ii).XYZ();
66 dist = Abs(a*xyz.X() + b*xyz.Y() + c*xyz.Z() + d);
67 if (dist > dfMaxDist)
68 dfMaxDist = dist;
69 }
70
71 return dfMaxDist;
72}
73
74//=======================================================================
75//function : BRepLib_FindSurface
76//purpose :
77//=======================================================================
78BRepLib_FindSurface::BRepLib_FindSurface()
79{
80}
81//=======================================================================
82//function : BRepLib_FindSurface
83//purpose :
84//=======================================================================
85BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape& S,
86 const Standard_Real Tol,
87 const Standard_Boolean OnlyPlane)
88{
89 Init(S,Tol,OnlyPlane);
90}
91//=======================================================================
92//function : Init
93//purpose :
94//=======================================================================
95void BRepLib_FindSurface::Init(const TopoDS_Shape& S,
96 const Standard_Real Tol,
97 const Standard_Boolean OnlyPlane)
98{
99 myTolerance = Tol;
100 myTolReached = 0.;
101 isExisted = Standard_False;
102 myLocation.Identity();
103 mySurface.Nullify();
104
105 // compute the tolerance
106 TopExp_Explorer ex;
107
108 for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
109 Standard_Real t = BRep_Tool::Tolerance(TopoDS::Edge(ex.Current()));
110 if (t > myTolerance) myTolerance = t;
111 }
112
113 // search an existing surface
114 ex.Init(S,TopAbs_EDGE);
115 if (!ex.More()) return; // no edges ....
116
117 TopoDS_Edge E = TopoDS::Edge(ex.Current());
118 Standard_Real f,l,ff,ll;
119 Handle(Geom2d_Curve) PC,PPC;
120 Handle(Geom_Surface) SS;
121 TopLoc_Location L;
122 Standard_Integer i = 0,j;
123
124 // iterate on the surfaces of the first edge
125 while ( Standard_True) {
126 i++;
127 BRep_Tool::CurveOnSurface(E,PC,mySurface,myLocation,f,l,i);
128 if (mySurface.IsNull()) {
129 break;
130 }
131 // check the other edges
132 for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
133 if (!E.IsSame(ex.Current())) {
134 j = 0;
135 while (Standard_True) {
136 j++;
137 BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()),
138 PPC,SS,L,ff,ll,j);
139 if (SS.IsNull()) {
140 break;
141 }
142 if (SS == mySurface) {
143 break;
144 }
145 SS.Nullify();
146 }
147
148 if (SS.IsNull()) {
149 mySurface.Nullify();
150 break;
151 }
152 }
153 }
154
155 // if OnlyPlane, eval if mySurface is a plane.
156 if ( OnlyPlane && !mySurface.IsNull() )
157 mySurface = Handle(Geom_Plane)::DownCast(mySurface);
158
159 if (!mySurface.IsNull()) break;
160 }
161
162 if (!mySurface.IsNull()) {
163 isExisted = Standard_True;
164 return;
165 }
166 //
167 // no existing surface, search a plane
168 // 07/02/02 akm vvv : (OCC157) changed algorithm
169 // 1. Collect the points along all edges of the shape
170 // For each point calculate the WEIGHT = sum of
171 // distances from neighboring points (_only_ same edge)
172 // 2. Minimizing the weighed sum of squared deviations
173 // compute coefficients of the sought plane.
174
175 TColgp_SequenceOfPnt aPoints;
176 TColStd_SequenceOfReal aWeight;
177
178 // ======================= Step #1
179 for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next())
180 {
181 BRepAdaptor_Curve c(TopoDS::Edge(ex.Current()));
182
183 Standard_Real dfUf = c.FirstParameter();
184 Standard_Real dfUl = c.LastParameter();
185 if (IsEqual(dfUf,dfUl)) {
186 // Degenerate
187 continue;
188 }
189 Standard_Integer iNbPoints=0;
190
191 // Add the points with weights to the sequences
192 switch (c.GetType())
193 {
194 case GeomAbs_BezierCurve:
195 {
196 // Put all poles for bezier
197 Handle(Geom_BezierCurve) GC = c.Bezier();
198 Standard_Integer iNbPol = GC->NbPoles();
199 if ( iNbPol < 2)
200 // Degenerate
201 continue;
202 else
203 {
204 Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
205 GC->Poles(aPoles->ChangeArray1());
206 gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
207 Standard_Real dfDistPrev = 0., dfDistNext;
208 for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
209 {
210 if (iPol<iNbPol)
211 {
212 aPoleNext = aPoles->Value(iPol+1);
213 dfDistNext = aPolePrev.Distance(aPoleNext);
214 }
215 else
216 dfDistNext = 0.;
217 aPoints.Append (aPolePrev);
218 aWeight.Append (dfDistPrev+dfDistNext);
219 dfDistPrev = dfDistNext;
220 aPolePrev = aPoleNext;
221 }
222 }
223 }
224 break;
225 case GeomAbs_BSplineCurve:
226 {
227 // Put all poles for bspline
228 Handle(Geom_BSplineCurve) GC = c.BSpline();
229 Standard_Integer iNbPol = GC->NbPoles();
230 if ( iNbPol < 2)
231 // Degenerate
232 continue;
233 else
234 {
235 Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
236 GC->Poles(aPoles->ChangeArray1());
237 gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
238 Standard_Real dfDistPrev = 0., dfDistNext;
239 for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
240 {
241 if (iPol<iNbPol)
242 {
243 aPoleNext = aPoles->Value(iPol+1);
244 dfDistNext = aPolePrev.Distance(aPoleNext);
245 }
246 else
247 dfDistNext = 0.;
248 aPoints.Append (aPolePrev);
249 aWeight.Append (dfDistPrev+dfDistNext);
250 dfDistPrev = dfDistNext;
251 aPolePrev = aPoleNext;
252 }
253 }
254 }
255 break;
256
257 case GeomAbs_Line:
258 case GeomAbs_Circle:
259 case GeomAbs_Ellipse:
260 case GeomAbs_Hyperbola:
261 case GeomAbs_Parabola:
262 if (c.GetType() == GeomAbs_Line)
263 // Two points on straight segment
264 iNbPoints=2;
265 else
266 // Four points on otheranalitical curves
267 iNbPoints=4;
268 default:
269 {
270 // Put some points on other curves
271 if (iNbPoints==0)
272 iNbPoints = 15 + c.NbIntervals(GeomAbs_C3);
273 Standard_Real dfDelta = (dfUl-dfUf)/(iNbPoints-1);
274 Standard_Integer iPoint;
275 Standard_Real dfU;
276 gp_Pnt aPointPrev = c.Value(dfUf), aPointNext;
277 Standard_Real dfDistPrev = 0., dfDistNext;
278 for (iPoint=1, dfU=dfUf+dfDelta;
279 iPoint<=iNbPoints;
280 iPoint++, dfU+=dfDelta)
281 {
282 if (iPoint<iNbPoints)
283 {
284 aPointNext = c.Value(dfU);
285 dfDistNext = aPointPrev.Distance(aPointNext);
286 }
287 else
288 dfDistNext = 0.;
289 aPoints.Append (aPointPrev);
290 aWeight.Append (dfDistPrev+dfDistNext);
291 dfDistPrev = dfDistNext;
292 aPointPrev = aPointNext;
293 }
294 } // default:
295 } // switch (c.GetType()) ...
296 } // for (ex.Init(S,TopAbs_EDGE); ex.More() && control; ex.Next()) ...
297
298 if (aPoints.Length() < 3) {
299 return;
300 }
301
302 // ======================= Step #2
303 myLocation.Identity();
304 Standard_Integer iPoint;
305 math_Matrix aMat (1,3,1,3, 0.);
306 math_Vector aVec (1,3, 0.);
307 // Find the barycenter and normalize weights
308 Standard_Real dfMaxWeight=0.;
309 gp_XYZ aBaryCenter(0.,0.,0.);
310 Standard_Real dfSumWeight=0.;
311 for (iPoint=1; iPoint<=aPoints.Length(); iPoint++) {
312 Standard_Real dfW = aWeight(iPoint);
313 aBaryCenter += dfW*aPoints(iPoint).XYZ();
314 dfSumWeight += dfW;
315 if (dfW > dfMaxWeight) {
316 dfMaxWeight = dfW;
317 }
318 }
319 aBaryCenter /= dfSumWeight;
320
321 // Fill the matrix and the right vector
322 for (iPoint=1; iPoint<=aPoints.Length(); iPoint++) {
323 gp_XYZ p=aPoints(iPoint).XYZ()-aBaryCenter;
324 Standard_Real w=aWeight(iPoint)/dfMaxWeight;
325 aMat(1,1)+=w*p.X()*p.X();
326 aMat(1,2)+=w*p.X()*p.Y();
327 aMat(1,3)+=w*p.X()*p.Z();
328 aMat(2,1)+=w*p.Y()*p.X();
329 aMat(2,2)+=w*p.Y()*p.Y();
330 aMat(2,3)+=w*p.Y()*p.Z();
331 aMat(3,1)+=w*p.Z()*p.X();
332 aMat(3,2)+=w*p.Z()*p.Y();
333 aMat(3,3)+=w*p.Z()*p.Z();
334 aVec(1) -= w*p.X();
335 aVec(2) -= w*p.Y();
336 aVec(3) -= w*p.Z();
337 }
338
339 // Solve the system of equations to get plane coefficients
340 math_Gauss aSolver(aMat);
341 Standard_Boolean isSolved = aSolver.IsDone();
342 //
343 // let us be more tolerant (occ415)
344 Standard_Real dfDist = RealLast();
345 Handle(Geom_Plane) aPlane;
346 //
347 if (isSolved) {
348 aSolver.Solve(aVec);
349 if (aVec.Norm2()<gp::Resolution()) {
350 isSolved = Standard_False;
351 }
352 }
353 //
354 if (isSolved) {
355 aPlane = new Geom_Plane(aBaryCenter,gp_Dir(aVec(1),aVec(2),aVec(3)));
356 dfDist = Controle (aPoints, aPlane);
357 }
358 //
359 if (!isSolved || myTolerance < dfDist) {
360 gp_Pnt aFirstPnt=aPoints(1);
361 for (iPoint=2; iPoint<=aPoints.Length(); iPoint++) {
362 gp_Vec aDir(aFirstPnt,aPoints(iPoint));
363 Standard_Real dfSide=aDir.Magnitude();
364 if (dfSide<myTolerance) {
365 continue; // degeneration
366 }
367 for (Standard_Integer iP1=iPoint+1; iP1<=aPoints.Length(); iP1++) {
368 gp_Vec aCross = gp_Vec(aFirstPnt,aPoints(iP1)) ^ aDir ;
369 if (aCross.Magnitude() > dfSide*myTolerance) {
370 Handle(Geom_Plane) aPlane2 = new Geom_Plane(aFirstPnt, aCross);
371 Standard_Real dfDist2 = Controle (aPoints, aPlane2);
372 if (dfDist2 < myTolerance) {
373 myTolReached = dfDist2;
374 mySurface = aPlane2;
375 return;
376 }
377 if (dfDist2 < dfDist) {
378 dfDist = dfDist2;
379 aPlane = aPlane2;
380 }
381 }
382 }
383 }
384 }
385 //
386 //XXf
387 //static Standard_Real weakness = 5.0;
388 Standard_Real weakness = 5.0;
389 //XXf
390 if(dfDist <= myTolerance || dfDist < myTolerance*weakness && Tol<0) {
391 //XXf
392 //myTolReached = dfDist;
393 //XXt
394 mySurface = aPlane;
395 }
396 //XXf
397 myTolReached = dfDist;
398 //XXt
399}
400//=======================================================================
401//function : Found
402//purpose :
403//=======================================================================
404Standard_Boolean BRepLib_FindSurface::Found() const
405{
406 return !mySurface.IsNull();
407}
408//=======================================================================
409//function : Surface
410//purpose :
411//=======================================================================
412Handle(Geom_Surface) BRepLib_FindSurface::Surface() const
413{
414 return mySurface;
415}
416//=======================================================================
417//function : Tolerance
418//purpose :
419//=======================================================================
420Standard_Real BRepLib_FindSurface::Tolerance() const
421{
422 return myTolerance;
423}
424//=======================================================================
425//function : ToleranceReached
426//purpose :
427//=======================================================================
428Standard_Real BRepLib_FindSurface::ToleranceReached() const
429{
430 return myTolReached;
431}
432//=======================================================================
433//function : Existed
434//purpose :
435//=======================================================================
436Standard_Boolean BRepLib_FindSurface::Existed() const
437{
438 return isExisted;
439}
440//=======================================================================
441//function : Location
442//purpose :
443//=======================================================================
444TopLoc_Location BRepLib_FindSurface::Location() const
445{
446 return myLocation;
447}
Open CASCADE Technology repository