0031007: Coding - eliminate warnings issued while compiling with -pedantic flag
[occt.git] / src / IntCurveSurface / IntCurveSurface_Polyhedron.gxx
1 // Created on: 1993-02-03
2 // Created by: Laurent BUCHARD
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 #include <gp.hxx>
18 #include <gp_Pnt.hxx>
19 #include <gp_Vec.hxx>
20 #include <gp_Dir.hxx>
21 #include <gp_Lin.hxx>
22 #include <TColgp_Array2OfPnt.hxx>
23 #include <TColStd_Array2OfReal.hxx>
24
25 #include <Bnd_Array1OfBox.hxx>
26
27 #include <Standard_ConstructionError.hxx>
28
29 #include <stdio.h>
30
31 #define CHECKBOUNDS 0
32
33
34 //-----------------------------------------------------
35 #define LONGUEUR_MINI_EDGE_TRIANGLE 1e-15
36
37
38 //=======================================================================
39 //function : IntCurveSurface_Polyhedron
40 //purpose  : 
41 //=======================================================================
42 IntCurveSurface_Polyhedron::IntCurveSurface_Polyhedron (const ThePSurface&     Surface,
43                                                         const Standard_Integer nbdU,
44                                                         const Standard_Integer nbdV,
45                                                         const Standard_Real    u1,
46                                                         const Standard_Real    v1,
47                                                         const Standard_Real    u2,
48                                                         const Standard_Real    v2)
49      : nbdeltaU((nbdU<3)? 3 : nbdU),
50        nbdeltaV((nbdV<3)? 3 : nbdV),
51        TheDeflection(Epsilon(100.)),
52        C_MyPnts(NULL),C_MyU(NULL),C_MyV(NULL),C_MyIsOnBounds(NULL)
53 {
54   Standard_Integer t = (nbdeltaU+1)*(nbdeltaV+1)+1; 
55   gp_Pnt           *CMyPnts = new gp_Pnt[t];         C_MyPnts = (void *)CMyPnts;
56   Standard_Real    *CMyU    = new Standard_Real[t];  C_MyU    = (void *)CMyU;  
57   Standard_Real    *CMyV    = new Standard_Real[t];  C_MyV    = (void *)CMyV;
58
59 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:46 2001 Begin
60   Standard_Boolean *CMyIsOnBounds = new Standard_Boolean[t];
61
62   C_MyIsOnBounds = (void *)CMyIsOnBounds;
63 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:47 2001 End
64   Init(Surface,u1,v1,u2,v2);
65 }
66
67 //=======================================================================
68 //function : IntCurveSurface_Polyhedron
69 //purpose  : 
70 //=======================================================================
71 IntCurveSurface_Polyhedron::IntCurveSurface_Polyhedron (const ThePSurface&     Surface,
72                                                         const TColStd_Array1OfReal&    Upars,
73                                                         const TColStd_Array1OfReal&    Vpars)
74     : nbdeltaU(Upars.Length()-1),
75       nbdeltaV(Vpars.Length()-1),
76       TheDeflection(Epsilon(100.)),
77       C_MyPnts(NULL),C_MyU(NULL),C_MyV(NULL),C_MyIsOnBounds(NULL)
78 {
79   Standard_Integer t = (nbdeltaU+1)*(nbdeltaV+1)+1; 
80   gp_Pnt           *CMyPnts = new gp_Pnt[t];         C_MyPnts = (void *)CMyPnts;
81   Standard_Real    *CMyU    = new Standard_Real[t];  C_MyU    = (void *)CMyU;  
82   Standard_Real    *CMyV    = new Standard_Real[t];  C_MyV    = (void *)CMyV;
83
84 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:46 2001 Begin
85   Standard_Boolean *CMyIsOnBounds = new Standard_Boolean[t];
86
87   C_MyIsOnBounds = (void *)CMyIsOnBounds;
88 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:47 2001 End
89   Init(Surface, Upars, Vpars);
90 }
91
92
93 void IntCurveSurface_Polyhedron::Destroy() { 
94   //-- printf("\n IntCurveSurface_Polyhedron::Destroy\n");
95   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;       if(C_MyPnts) delete [] CMyPnts;
96   Standard_Real *CMyU = (Standard_Real *)C_MyU;   if(C_MyU)    delete [] CMyU;
97   Standard_Real *CMyV = (Standard_Real *)C_MyV;   if(C_MyV)    delete [] CMyV;
98
99 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:46 2001 Begin
100   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
101
102   if(C_MyIsOnBounds)    delete [] CMyIsOnBounds;
103 //  Modified by Sergey KHROMOV - Fri Dec  7 12:03:47 2001 End
104
105   C_MyPnts=C_MyU=C_MyV=C_MyIsOnBounds=NULL;
106 }
107
108 //=======================================================================
109 //function : Init
110 //purpose  : 
111 //=======================================================================
112 void IntCurveSurface_Polyhedron::Init(const ThePSurface&     Surface,
113                                       const Standard_Real U0,
114                                       const Standard_Real V0,
115                                       const Standard_Real U1,
116                                       const Standard_Real V1) { 
117   //#define DEBUGDUMP
118   Standard_Integer i1,i2;
119   Standard_Real    U,V;
120   Standard_Real    U1mU0sNbdeltaU = (U1-U0)/(Standard_Real)nbdeltaU;
121   Standard_Real    V1mV0sNbdeltaV = (V1-V0)/(Standard_Real)nbdeltaV;
122   gp_Pnt TP;
123   Standard_Integer Index=1;
124   //-- --------------------------------------------------
125   //-- Index varie de 1 -> (nbdu+1)*(nbdv+1)
126   //-- V est la colonne
127   //-- U est la ligne 
128   //-- --------------------------------------------------
129   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
130   Standard_Real *CMyU = (Standard_Real *)C_MyU;
131   Standard_Real *CMyV = (Standard_Real *)C_MyV;
132   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
133
134   for (i1 = 0, U = U0; i1 <= nbdeltaU; i1++, U+= U1mU0sNbdeltaU) {
135     for (i2 = 0, V = V0; i2 <= nbdeltaV; i2++, V+= V1mV0sNbdeltaV ) {
136       ThePSurfaceTool::D0(Surface,U,V,TP);
137       //-- Point(TP,i1, i2,U,V);
138       CMyPnts[Index] = TP;
139       CMyU[Index]    = U;
140       CMyV[Index]    = V;
141 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:51 2001
142       CMyIsOnBounds[Index] = (i1 == 0 || i1 == nbdeltaU ||
143                               i2 == 0 || i2 == nbdeltaV);
144 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:52 2001
145       TheBnd.Add(TP);
146       Index++;
147     }
148   }
149   //-- Calcul de la deflection Triangle <-> Point milieu
150   Standard_Real tol=0.0; Standard_Integer nbtriangles = NbTriangles();
151   for (i1=1; i1<=nbtriangles; i1++) {
152     Standard_Real tol1 = DeflectionOnTriangle(Surface,i1);
153     if(tol1>tol) tol=tol1;
154   }
155   //-- Calcul de la deflection Bord <-> Point milieu
156
157
158   DeflectionOverEstimation(tol*1.2);
159   FillBounding();
160
161 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:33 2001 Begin
162   Standard_Real aDeflection;
163
164   TheBorderDeflection = RealFirst();
165
166 // Compute the deflection on the lower bound (U-isoline) of the surface.
167   aDeflection = ComputeBorderDeflection(Surface, U0, V0, V1, Standard_True);
168
169   if (aDeflection > TheBorderDeflection)
170     TheBorderDeflection = aDeflection;
171
172 // Compute the deflection on the upper bound (U-isoline) of the surface.
173   aDeflection = ComputeBorderDeflection(Surface, U1, V0, V1, Standard_True);
174
175   if (aDeflection > TheBorderDeflection)
176     TheBorderDeflection = aDeflection;
177
178 // Compute the deflection on the lower bound (V-isoline) of the surface.
179   aDeflection = ComputeBorderDeflection(Surface, V0, U0, U1, Standard_False);
180
181   if (aDeflection > TheBorderDeflection)
182     TheBorderDeflection = aDeflection;
183
184 // Compute the deflection on the upper bound (V-isoline) of the surface.
185   aDeflection = ComputeBorderDeflection(Surface, V1, U0, U1, Standard_False);
186
187   if (aDeflection > TheBorderDeflection)
188     TheBorderDeflection = aDeflection;
189
190 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:34 2001 End
191
192   #ifdef OCCT_DEBUG_DUMP 
193     Dump();
194   #endif
195 }
196 //=======================================================================
197 //function : Init
198 //purpose  : 
199 //=======================================================================
200 void IntCurveSurface_Polyhedron::Init(const ThePSurface&     Surface,
201                                       const TColStd_Array1OfReal& Upars,
202                                       const TColStd_Array1OfReal& Vpars) { 
203   //#define DEBUGDUMP
204   Standard_Integer i1,i2;
205   Standard_Real    U,V;
206   gp_Pnt TP;
207   Standard_Integer Index=1;
208   //-- --------------------------------------------------
209   //-- Index varie de 1 -> (nbdu+1)*(nbdv+1)
210   //-- V est la colonne
211   //-- U est la ligne 
212   //-- --------------------------------------------------
213   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
214   Standard_Real *CMyU = (Standard_Real *)C_MyU;
215   Standard_Real *CMyV = (Standard_Real *)C_MyV;
216   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
217   Standard_Integer i0 = Upars.Lower(), j0 = Vpars.Lower();
218
219   for (i1 = 0; i1 <= nbdeltaU; i1++) {
220     U = Upars(i1+i0);
221     for (i2 = 0; i2 <= nbdeltaV; i2++) {
222       V = Vpars(i2+j0);
223       ThePSurfaceTool::D0(Surface,U,V,TP);
224       //-- Point(TP,i1, i2,U,V);
225       CMyPnts[Index] = TP;
226       CMyU[Index]    = U;
227       CMyV[Index]    = V;
228 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:51 2001
229       CMyIsOnBounds[Index] = (i1 == 0 || i1 == nbdeltaU ||
230                               i2 == 0 || i2 == nbdeltaV);
231 //  Modified by Sergey KHROMOV - Fri Dec  7 12:07:52 2001
232       TheBnd.Add(TP);
233       Index++;
234     }
235   }
236   //-- Calcul de la deflection Triangle <-> Point milieu
237   Standard_Real tol=0.0; Standard_Integer nbtriangles = NbTriangles();
238   for (i1=1; i1<=nbtriangles; i1++) {
239     Standard_Real tol1 = DeflectionOnTriangle(Surface,i1);
240     if(tol1>tol) tol=tol1;
241   }
242   //-- Calcul de la deflection Bord <-> Point milieu
243
244
245   DeflectionOverEstimation(tol*1.2);
246   FillBounding();
247
248 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:33 2001 Begin
249   Standard_Real aDeflection;
250
251   TheBorderDeflection = RealFirst();
252   Standard_Real U0 = Upars(i0);
253   Standard_Real V0 = Vpars(j0);
254   Standard_Real U1 = Upars(Upars.Upper());
255   Standard_Real V1 = Vpars(Vpars.Upper());
256
257 // Compute the deflection on the lower bound (U-isoline) of the surface.
258   aDeflection = ComputeBorderDeflection(Surface, U0, V0, V1, Standard_True);
259
260   if (aDeflection > TheBorderDeflection)
261     TheBorderDeflection = aDeflection;
262
263 // Compute the deflection on the upper bound (U-isoline) of the surface.
264   aDeflection = ComputeBorderDeflection(Surface, U1, V0, V1, Standard_True);
265
266   if (aDeflection > TheBorderDeflection)
267     TheBorderDeflection = aDeflection;
268
269 // Compute the deflection on the lower bound (V-isoline) of the surface.
270   aDeflection = ComputeBorderDeflection(Surface, V0, U0, U1, Standard_False);
271
272   if (aDeflection > TheBorderDeflection)
273     TheBorderDeflection = aDeflection;
274
275 // Compute the deflection on the upper bound (V-isoline) of the surface.
276   aDeflection = ComputeBorderDeflection(Surface, V1, U0, U1, Standard_False);
277
278   if (aDeflection > TheBorderDeflection)
279     TheBorderDeflection = aDeflection;
280
281 //  Modified by Sergey KHROMOV - Fri Dec  7 11:23:34 2001 End
282
283   #ifdef OCCT_DEBUG_DUMP 
284     Dump();
285   #endif
286 }
287 //=======================================================================
288 //function : DeflectionOnTriangle
289 //purpose  : 
290 //=======================================================================
291 Standard_Real IntCurveSurface_Polyhedron::DeflectionOnTriangle (const ThePSurface& Surface,
292                                                                 const Standard_Integer Triang) const 
293 {
294   Standard_Integer i1,i2,i3;    
295   Triangle(Triang,i1,i2,i3);
296   //-- Calcul de l equation du plan
297   Standard_Real u1,v1,u2,v2,u3,v3;
298   gp_Pnt P1,P2,P3;
299   P1 = Point(i1,u1,v1);
300   P2 = Point(i2,u2,v2);
301   P3 = Point(i3,u3,v3);
302   if(P1.SquareDistance(P2)<=LONGUEUR_MINI_EDGE_TRIANGLE) return(0);
303   if(P1.SquareDistance(P3)<=LONGUEUR_MINI_EDGE_TRIANGLE) return(0);
304   if(P2.SquareDistance(P3)<=LONGUEUR_MINI_EDGE_TRIANGLE) return(0);
305   gp_XYZ XYZ1=P2.XYZ()-P1.XYZ();
306   gp_XYZ XYZ2=P3.XYZ()-P2.XYZ();
307   gp_XYZ XYZ3=P1.XYZ()-P3.XYZ();
308   gp_Vec NormalVector((XYZ1^XYZ2)+(XYZ2^XYZ3)+(XYZ3^XYZ1));
309   Standard_Real aNormLen = NormalVector.Magnitude();
310   if (aNormLen < gp::Resolution()) {
311     return 0.;
312   }
313   //
314   NormalVector.Divide(aNormLen);
315   //-- Standard_Real PolarDistance = NormalVector * P1.XYZ();
316   //-- Calcul du point u,v  au centre du triangle
317   Standard_Real u = (u1+u2+u3)/3.0;
318   Standard_Real v = (v1+v2+v3)/3.0;
319   gp_Pnt P =  ThePSurfaceTool::Value(Surface,u,v);
320   gp_Vec P1P(P1,P);
321   return(Abs(P1P.Dot(NormalVector)));
322 }
323 //=======================================================================
324 //function : Parameters
325 //purpose  : 
326 //=======================================================================
327 void IntCurveSurface_Polyhedron::Parameters( const Standard_Integer Index
328                                             ,Standard_Real &U
329                                             ,Standard_Real &V) const 
330 {
331 #if CHECKBOUNDS 
332   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) { 
333     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
334   }
335 #endif
336   Standard_Real *CMyU = (Standard_Real *)C_MyU;
337   U = CMyU[Index];
338   Standard_Real *CMyV = (Standard_Real *)C_MyV;
339   V = CMyV[Index];
340 }
341 //=======================================================================
342 //function : DeflectionOverEstimation
343 //purpose  : Set
344 //=======================================================================
345 void IntCurveSurface_Polyhedron::DeflectionOverEstimation(const Standard_Real flec)
346 {
347   if(flec<0.0001) {  
348     TheDeflection=0.0001;
349     TheBnd.Enlarge(0.0001);
350   }
351   else { 
352     TheDeflection=flec;
353     TheBnd.Enlarge(flec);
354   }    
355 }
356 //=======================================================================
357 //function : DeflectionOverEstimation
358 //purpose  : Get
359 //=======================================================================
360 Standard_Real IntCurveSurface_Polyhedron::DeflectionOverEstimation() const
361 {
362   return TheDeflection;
363 }
364 //=======================================================================
365 //function : Bounding
366 //purpose  : 
367 //=======================================================================
368 const Bnd_Box& IntCurveSurface_Polyhedron::Bounding() const
369 {
370   return TheBnd;
371 }
372 //=======================================================================
373 //function : FillBounding
374 //purpose  :  
375 //=======================================================================
376 void IntCurveSurface_Polyhedron::FillBounding()
377 {
378   TheComponentsBnd=new Bnd_HArray1OfBox(1, NbTriangles());
379   Bnd_Box Boite;
380   Standard_Integer np1, np2, np3;
381   Standard_Integer nbtriangles = NbTriangles();
382   for (Standard_Integer iTri=1; iTri<=nbtriangles; iTri++) {
383     Triangle(iTri, np1, np2, np3);
384     gp_Pnt p1(Point(np1));
385     gp_Pnt p2(Point(np2));
386     gp_Pnt p3(Point(np3));
387     Boite.SetVoid();
388     if(p1.SquareDistance(p2)>LONGUEUR_MINI_EDGE_TRIANGLE) {
389       if(p1.SquareDistance(p3)>LONGUEUR_MINI_EDGE_TRIANGLE) {
390         if(p2.SquareDistance(p3)>LONGUEUR_MINI_EDGE_TRIANGLE) {
391           Boite.Add(p1);
392           Boite.Add(p2);
393           Boite.Add(p3);
394           Boite.Enlarge(TheDeflection);
395         }  
396       }
397     }
398     Boite.Enlarge(TheDeflection);
399     TheComponentsBnd->SetValue(iTri,Boite); 
400   }
401 }
402 //=======================================================================
403 //function : ComponentsBounding
404 //purpose  : 
405 //=======================================================================
406 const Handle(Bnd_HArray1OfBox)& 
407       IntCurveSurface_Polyhedron::ComponentsBounding() const
408 {
409  return TheComponentsBnd;
410 }
411 //=======================================================================
412 //function : NbTriangles
413 //purpose  : 
414 //=======================================================================
415 Standard_Integer IntCurveSurface_Polyhedron::NbTriangles  () const
416 {
417   return nbdeltaU*nbdeltaV*2;
418 }
419 //=======================================================================
420 //function : NbPoints
421 //purpose  : 
422 //=======================================================================
423 Standard_Integer IntCurveSurface_Polyhedron::NbPoints () const
424 {
425   return (nbdeltaU+1)*(nbdeltaV+1);
426 }
427 //=======================================================================
428 //function : TriConnex
429 //purpose  : 
430 //=======================================================================
431 Standard_Integer IntCurveSurface_Polyhedron::TriConnex
432   (const Standard_Integer Triang,
433    const Standard_Integer Pivot,
434    const Standard_Integer Pedge,
435    Standard_Integer&      TriCon,
436    Standard_Integer&      OtherP) const
437 {
438   Standard_Integer Pivotm1    = Pivot-1;
439   Standard_Integer nbdeltaVp1 = nbdeltaV+1;
440   Standard_Integer nbdeltaVm2 = nbdeltaV + nbdeltaV;
441
442 // Pivot position in the MaTriangle :
443   Standard_Integer ligP = Pivotm1/nbdeltaVp1;
444   Standard_Integer colP = Pivotm1 - ligP * nbdeltaVp1;
445
446 // Point sur Edge position in the MaTriangle and edge typ :
447   Standard_Integer ligE =0, colE =0, typE =0;
448   if (Pedge!=0) {
449     ligE= (Pedge-1)/nbdeltaVp1;
450     colE= (Pedge-1) - (ligE * nbdeltaVp1);
451   // Horizontal
452     if      (ligP==ligE) typE=1;
453   // Vertical
454     else if (colP==colE) typE=2;
455   // Oblique
456     else                 typE=3;
457   }
458   else {
459     typE=0;
460   }
461
462 // Triangle position General case :
463
464   Standard_Integer linT =0, colT =0;
465   Standard_Integer linO =0, colO =0;
466   Standard_Integer t =0, tt =0;
467
468   if (Triang!=0) {
469     t = (Triang-1)/(nbdeltaVm2);
470     tt= (Triang-1)-t*nbdeltaVm2;
471     linT= 1+t;
472     colT= 1+tt;
473     if (typE==0) {
474       if (ligP==linT) {
475         ligE=ligP-1;
476         colE=colP-1;
477         typE=3;
478       }
479       else {
480         if (colT==ligP+ligP) {
481           ligE=ligP;
482           colE=colP-1;
483           typE=1;
484         }
485         else {
486           ligE=ligP+1;
487           colE=colP+1;
488           typE=3;
489         }
490       }
491     }
492     switch (typE) {
493     case 1:  // Horizontal
494       if (linT==ligP) {
495         linT++;
496         linO=ligP+1;
497         colO=(colP>colE)? colP : colE; //--colO=Max(colP, colE);
498       }
499       else {
500         linT--;
501         linO=ligP-1;
502         colO=(colP<colE)? colP : colE;  //--colO=Min(colP, colE);
503       }
504       break;
505     case 2:  // Vertical
506       if (colT==(colP+colP)) {
507         colT++;
508         linO=(ligP>ligE)? ligP : ligE;  //--linO=Max(ligP, ligE);
509         colO=colP+1;
510       }
511       else {
512         colT--;
513         linO=(ligP<ligE)? ligP : ligE;  //--linO=Min(ligP, ligE);
514         colO=colP-1;
515       }
516       break;
517     case 3:  // Oblique
518       if ((colT&1)==0) {
519         colT--;
520         linO=(ligP>ligE)? ligP : ligE;  //--linO=Max(ligP, ligE);
521         colO=(colP<colE)? colP : colE;  //--colO=Min(colP, colE);
522       }
523       else {
524         colT++;
525         linO=(ligP<ligE)? ligP : ligE;  //--linO=Min(ligP, ligE);
526         colO=(colP>colE)? colP : colE;  //--colO=Max(colP, colE);
527       }
528       break;
529     }
530   }
531   else {
532     // Unknown Triangle position :
533     if (Pedge==0) {
534       // Unknown edge :
535       linT=(1>ligP)? 1 : ligP;      //--linT=Max(1, ligP);
536       colT=(1>(colP+colP))? 1 : (colP+colP);       //--colT=Max(1, colP+colP);
537       if (ligP==0) linO=ligP+1;
538       else         linO=ligP-1;
539       colO=colP;
540     }
541     else {
542       // Known edge We take the left or down connectivity :
543       switch (typE) {
544       case 1:  // Horizontal
545         linT=ligP+1;
546         colT=(colP>colE)? colP : colE;  //--colT=Max(colP,colE);
547         colT+=colT;
548         linO=ligP+1;
549         colO=(colP>colE)? colP : colE;  //--colO=Max(colP,colE);
550         break;
551       case 2:  // Vertical
552         linT=(ligP>ligE)? ligP : ligE;  //--linT=Max(ligP, ligE);
553         colT=colP+colP;
554         linO=(ligP<ligE)? ligP : ligE;  //--linO=Min(ligP, ligE);
555         colO=colP-1;
556         break;
557       case 3:  // Oblique
558         linT=(ligP>ligE)? ligP : ligE;  //--linT=Max(ligP, ligE);
559         colT=colP+colE;
560         linO=(ligP>ligE)? ligP : ligE;  //--linO=Max(ligP, ligE);
561         colO=(colP<colE)? colP : colE;  //--colO=Min(colP, colE);
562         break;
563       }
564     }
565   }
566
567   TriCon=(linT-1)*nbdeltaVm2 + colT;
568
569   if (linT<1) {
570     linO=0;
571     colO=colP+colP-colE;
572     if (colO<0) {colO=0;linO=1;}
573     else if (colO>nbdeltaV) {colO=nbdeltaV;linO=1;}
574     TriCon=0;
575   }
576   else if (linT>nbdeltaU) {
577     linO=nbdeltaU;
578     colO=colP+colP-colE;
579     if (colO<0) {colO=0;linO=nbdeltaU-1;}
580     else if (colO>nbdeltaV) {colO=nbdeltaV;linO=nbdeltaU-1;}
581     TriCon=0;
582   }
583
584   if (colT<1) {
585     colO=0;
586     linO=ligP+ligP-ligE;
587     if (linO<0) {linO=0;colO=1;}
588     else if (linO>nbdeltaU) {linO=nbdeltaU;colO=1;}
589     TriCon=0;
590   }
591   else if (colT>nbdeltaV) {
592     colO=nbdeltaV;
593     linO=ligP+ligP-ligE;
594     if (linO<0) {linO=0;colO=nbdeltaV-1;}
595     else if (linO>nbdeltaU) {linO=nbdeltaU;colO=nbdeltaV-1;}
596     TriCon=0;
597   }
598
599   OtherP=linO*nbdeltaVp1 + colO+1;
600
601   return TriCon;
602 }
603
604
605 //=======================================================================
606 //function : PlaneEquation
607 //purpose  : 
608 //=======================================================================
609 void IntCurveSurface_Polyhedron::PlaneEquation (const Standard_Integer Triang,
610                                         gp_XYZ&        NormalVector,
611                                         Standard_Real& PolarDistance)  const
612 {
613   Standard_Integer i1,i2,i3;
614   Triangle(Triang,i1,i2,i3);
615
616   //--  gp_XYZ v1=Point(i2).XYZ()-Point(i1).XYZ();
617   //--  gp_XYZ v2=Point(i3).XYZ()-Point(i2).XYZ();
618   //--  gp_XYZ v3=Point(i1).XYZ()-Point(i3).XYZ();
619
620   gp_XYZ Pointi1(Point(i1).XYZ());
621   gp_XYZ Pointi2(Point(i2).XYZ());
622   gp_XYZ Pointi3(Point(i3).XYZ());
623   
624
625   gp_XYZ v1= Pointi2 - Pointi1;
626   gp_XYZ v2= Pointi3 - Pointi2;
627   gp_XYZ v3= Pointi1 - Pointi3;
628
629   if(v1.SquareModulus()<=LONGUEUR_MINI_EDGE_TRIANGLE) { NormalVector.SetCoord(1.0,0.0,0.0); return; } 
630   if(v2.SquareModulus()<=LONGUEUR_MINI_EDGE_TRIANGLE) { NormalVector.SetCoord(1.0,0.0,0.0); return; } 
631   if(v3.SquareModulus()<=LONGUEUR_MINI_EDGE_TRIANGLE) { NormalVector.SetCoord(1.0,0.0,0.0); return; } 
632
633   NormalVector= (v1^v2)+(v2^v3)+(v3^v1);
634   Standard_Real aNormLen = NormalVector.Modulus();
635   if (aNormLen < gp::Resolution()) {
636     PolarDistance = 0.;
637   }
638   else {
639     NormalVector.Divide(aNormLen);
640     PolarDistance = NormalVector * Point(i1).XYZ();
641   }
642 }
643 //=======================================================================
644 //function : Contain
645 //purpose  : 
646 //=======================================================================
647 Standard_Boolean IntCurveSurface_Polyhedron::Contain (const Standard_Integer Triang,
648                                                       const gp_Pnt& ThePnt) const
649 {
650   Standard_Integer i1,i2,i3;
651   Triangle(Triang,i1,i2,i3);
652   gp_XYZ Pointi1(Point(i1).XYZ());
653   gp_XYZ Pointi2(Point(i2).XYZ());
654   gp_XYZ Pointi3(Point(i3).XYZ());
655   
656   gp_XYZ v1=(Pointi2-Pointi1)^(ThePnt.XYZ()-Pointi1);
657   gp_XYZ v2=(Pointi3-Pointi2)^(ThePnt.XYZ()-Pointi2);
658   gp_XYZ v3=(Pointi1-Pointi3)^(ThePnt.XYZ()-Pointi3);
659   if (v1*v2 >= 0. && v2*v3 >= 0. && v3*v1>=0.) 
660     return Standard_True;
661   else 
662     return Standard_False;
663 }
664 //=======================================================================
665 //function : Dump
666 //purpose  : 
667 //=======================================================================
668 void IntCurveSurface_Polyhedron::Dump() const
669 {
670
671 }
672 //=======================================================================
673 //function : Size
674 //purpose  : 
675 //=======================================================================
676 void IntCurveSurface_Polyhedron::Size (Standard_Integer& nbdu,
677                               Standard_Integer& nbdv) const
678 {
679   nbdu=nbdeltaU;
680   nbdv=nbdeltaV;
681 }
682 //=======================================================================
683 //function : Triangle
684 //purpose  : 
685 //=======================================================================
686 void IntCurveSurface_Polyhedron::Triangle (const Standard_Integer Index,
687                                   Standard_Integer& P1,
688                                   Standard_Integer& P2,
689                                   Standard_Integer& P3) const
690 {
691   Standard_Integer line=1+((Index-1)/(nbdeltaV*2));
692   Standard_Integer colon=1+((Index-1)%(nbdeltaV*2));
693   Standard_Integer colpnt=(colon+1)/2;
694
695 // General formula = (line-1)*(nbdeltaV+1)+colpnt
696   
697 //  Position of P1 = MesXYZ(line,colpnt);
698   P1= (line-1)*(nbdeltaV+1) + colpnt;
699
700 //  Position of P2= MesXYZ(line+1,colpnt+((colon-1)%2));
701   P2= line*(nbdeltaV+1) + colpnt+((colon-1)%2);
702
703 //  Position of P3= MesXYZ(line+(colon%2),colpnt+1);
704   P3= (line-1+(colon%2))*(nbdeltaV+1) + colpnt + 1;
705 }
706 //=======================================================================
707 //function : Point
708 //=======================================================================
709 const gp_Pnt& IntCurveSurface_Polyhedron::Point(const Standard_Integer Index
710                                                 ,Standard_Real& U
711                                                 ,Standard_Real& V) const 
712 {
713 #if CHECKBOUNDS 
714   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) { 
715     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
716   }
717 #endif
718   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
719   Standard_Real *CMyU = (Standard_Real *)C_MyU;
720   Standard_Real *CMyV = (Standard_Real *)C_MyV;
721   U=CMyU[Index];
722   V=CMyV[Index];
723   return CMyPnts[Index];
724 }
725 //=======================================================================
726 //function : Point
727 //=======================================================================
728 const gp_Pnt& IntCurveSurface_Polyhedron::Point(const Standard_Integer Index) const {
729 #if CHECKBOUNDS 
730   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) { 
731     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
732   }
733 #endif
734   
735   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
736   return CMyPnts[Index];
737 }
738 //=======================================================================
739 //function : Point
740 //=======================================================================
741 //void IntCurveSurface_Polyhedron::Point (const gp_Pnt& p, 
742 //                                      const Standard_Integer lig,
743 //                                      const Standard_Integer col,
744 //                                      const Standard_Real u,
745 //                                      const Standard_Real v) 
746 void IntCurveSurface_Polyhedron::Point (const gp_Pnt& , 
747                                         const Standard_Integer ,
748                                         const Standard_Integer ,
749                                         const Standard_Real ,
750                                         const Standard_Real ) 
751 {
752   printf("\n IntCurveSurface_Polyhedron::Point : Ne dois pas etre appelle\n");
753 }
754 //=======================================================================
755 //function : Point
756 //=======================================================================
757 void IntCurveSurface_Polyhedron::Point (const Standard_Integer Index,gp_Pnt& P) const
758 {
759 #if CHECKBOUNDS 
760   if(Index<0 || Index>((nbdeltaU+1)*(nbdeltaV+1))) { 
761     printf("\n Erreur IntCurveSurface_Polyhedron::Parameters\n");
762   }
763 #endif
764   
765   gp_Pnt *CMyPnts     = (gp_Pnt *)C_MyPnts;
766   P = CMyPnts[Index];
767 }
768
769 //  Modified by Sergey KHROMOV - Fri Dec  7 10:12:47 2001 Begin
770
771 //=======================================================================
772 //function : IsOnBound
773 //purpose  : This method returns true if the edge based on points with 
774 //           indices Index1 and Index2 represents a boundary edge.
775 //=======================================================================
776
777 Standard_Boolean IntCurveSurface_Polyhedron::IsOnBound
778                                  (const Standard_Integer Index1,
779                                   const Standard_Integer Index2) const
780 {
781 #if CHECKBOUNDS 
782   if(Index1<0 || Index1>((nbdeltaU+1)*(nbdeltaV+1))) { 
783     printf("\n Erreur IntCurveSurface_Polyhedron::IsOnBound\n");
784   }
785   if(Index2<0 || Index2>((nbdeltaU+1)*(nbdeltaV+1))) { 
786     printf("\n Erreur IntCurveSurface_Polyhedron::IsOnBound\n");
787   }
788 #endif
789   Standard_Boolean *CMyIsOnBounds = (Standard_Boolean *)C_MyIsOnBounds;
790   Standard_Integer  aDiff         = Abs(Index1 - Index2);
791   Standard_Integer  i;
792
793 // Check if points are neighbour ones.
794   if (aDiff != 1 && aDiff != nbdeltaV + 1)
795     return Standard_False;
796
797   for (i = 0; i <= nbdeltaU; i++) {
798     if ((Index1 == 1 + i*(nbdeltaV + 1)) && (Index2 == Index1 - 1))
799       return Standard_False;
800
801     if ((Index1 == (1 + i)*(nbdeltaV + 1)) && (Index2 == Index1 + 1))
802       return Standard_False;
803   }
804
805   return (CMyIsOnBounds[Index1] && CMyIsOnBounds[Index2]);
806 }
807
808 //=======================================================================
809 //function : ComputeBorderDeflection
810 //purpose  : This method computes and returns a deflection of isoline 
811 //           of given parameter on Surface.
812 //=======================================================================
813
814 Standard_Real IntCurveSurface_Polyhedron::ComputeBorderDeflection
815                               (const ThePSurface      &Surface,
816                                const Standard_Real     Parameter,
817                                const Standard_Real     PMin,
818                                const Standard_Real     PMax,
819                                const Standard_Boolean  isUIso) const
820 {
821   Standard_Integer aNbSamples;
822   Standard_Integer i;
823   
824   if (isUIso)
825     aNbSamples = nbdeltaV;
826   else
827     aNbSamples = nbdeltaU;
828
829   Standard_Real aDelta      = (PMax - PMin)/aNbSamples;
830   Standard_Real aPar        = PMin;
831   Standard_Real aDeflection = RealFirst();
832   gp_XYZ        aP1;
833   gp_XYZ        aP2;
834   gp_XYZ        aPMid;
835   gp_XYZ        aPParMid;
836
837   for (i = 0; i <= aNbSamples; i++, aPar+= aDelta) {
838     if (isUIso) {
839       aP1      = ThePSurfaceTool::Value(Surface, Parameter, aPar).XYZ();
840       aP2      = ThePSurfaceTool::Value(Surface, Parameter, aPar + aDelta).XYZ();
841       aPParMid = ThePSurfaceTool::Value(Surface, Parameter, aPar + aDelta/2.).XYZ();
842     } else {
843       aP1      = ThePSurfaceTool::Value(Surface, aPar,             Parameter).XYZ();
844       aP2      = ThePSurfaceTool::Value(Surface, aPar + aDelta,    Parameter).XYZ();
845       aPParMid = ThePSurfaceTool::Value(Surface, aPar + aDelta/2., Parameter).XYZ();
846     }
847     aPMid = (aP2 + aP1)/2.;
848
849     Standard_Real aDist = (aPMid - aPParMid).Modulus();
850
851     if (aDist > aDeflection)
852       aDeflection = aDist;
853   }
854
855   return aDeflection;
856 }
857
858 //  Modified by Sergey KHROMOV - Fri Dec  7 11:21:52 2001 End