summaryrefslogtreecommitdiff
path: root/inc/Extrema_GenExtCC.gxx
blob: ed1f3a4f165df9524769284272ac1e7fbf58a417 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
// File:	Extrema_GenExtCC.gxx
// Created:	Tue Jul 18 17:42:19 1995
// Author:	Modelistation
//		<model@metrox>

#include <StdFail_NotDone.hxx>
#include <math_FunctionSetRoot.hxx>
#include <math_NewtonFunctionSetRoot.hxx>
#include <TColStd_Array2OfInteger.hxx>
#include <TColStd_Array2OfReal.hxx>
#include <Standard_NullObject.hxx>
#include <Standard_OutOfRange.hxx>
#include <Precision.hxx>

Extrema_GenExtCC::Extrema_GenExtCC () : myDone (Standard_False)
{
}

Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1,
			      const Curve2& C2,
			      const Standard_Integer NbU, 
			      const Standard_Integer NbV,
			      const Standard_Real TolU, 
			      const Standard_Real TolV) : myF (C1,C2, Min(TolU, TolV)), myDone (Standard_False)
{
  SetCurveCache (1, new Cache (C1, C1.FirstParameter(), C1.LastParameter(), NbU, Standard_True));
  SetCurveCache (2, new Cache (C2, C2.FirstParameter(), C2.LastParameter(), NbV, Standard_True));
  Perform();
}


Extrema_GenExtCC::Extrema_GenExtCC (const Curve1& C1,
			      const Curve2& C2,
			      const Standard_Real Uinf,
			      const Standard_Real Usup,
			      const Standard_Real Vinf,
			      const Standard_Real Vsup,
			      const Standard_Integer NbU, 
			      const Standard_Integer NbV,
			      const Standard_Real TolU, 
			      const Standard_Real TolV) : myF (C1,C2), myDone (Standard_False)
{
  SetCurveCache (1, new Cache (C1, Uinf, Usup, NbU, Standard_True));
  SetCurveCache (2, new Cache (C2, Vinf, Vsup, NbV, Standard_True));
  Perform();
}

void Extrema_GenExtCC::SetCurveCache (const Standard_Integer theRank,
                                      const Handle(Cache)& theCache)
{
  Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_GenExtCC::SetCurveCache()")
  myF.SetCurve (theRank, *(Curve1*)theCache->CurvePtr());
  Standard_Integer anInd = theRank - 1;
  myCache[anInd] = theCache;
}

void Extrema_GenExtCC::SetTolerance (const Standard_Real Tol)
{
  myF.SetTolerance (Tol);
}


//=============================================================================
void Extrema_GenExtCC::Perform ()
/*-----------------------------------------------------------------------------
Fonction:
   Recherche de toutes les distances extremales entre les courbes C1 et C2.
  a partir de 2 echantillonnages (NbU,NbV).

Methode:
   L'algorithme part de l'hypothese que les echantillonnages sont suffisamment
  fins pour que, s'il existe N distances extremales entre les 2 courbes,
  alors il existe aussi N extrema entre les 2 ensembles de points.
  Ainsi, l'algorithme consiste a partir des extrema des echantillons
  pour trouver les extrema des courbes.
   Les extrema sont calcules par l'algorithme math_FunctionSetRoot avec les
  arguments suivants:
  - myF: Extrema_FuncExtCC cree a partir de C1 et C2,
  - UV: math_Vector dont les composantes sont les parametres des points de
    l'extremum sur les ensembles de points,
  - Tol: Min(TolU,TolV), (Prov.:math_FunctionSetRoot n'autorise pas un vecteur)
  - UVinf: math_Vector dont les composantes sont les bornes inferieures de u et
    v,
  - UVsup: math_Vector dont les composantes sont les bornes superieures de u et
    v.

Traitement:
  a- Constitution du tableau des square distances (TbDist2(0,NbU+1,0,NbV+1)):
      Le tableau est volontairement etendu; les lignes 0 et NbU+1 et les
     colonnes 0 et NbV+1 seront initialisees a RealFirst() ou RealLast()
     pour simplifier les tests effectues dans l'etape b
     (on n'a pas besoin de tester si le point est sur une extremite).
  b- Calcul des extrema:
      On recherche d'abord les minima et ensuite les maxima. Ces 2 traitements
     se passent de facon similaire:
  b.a- Initialisations:
      - des 'bords' du tableau TbDist2 (a RealLast() dans le cas des minima
        et a RealLast() dans le cas des maxima),
      - du tableau TbSel(0,NbU+1,0,NbV+1) de selection des points pour un
        calcul d'extremum local (a 0). Lorsqu'un couple de points sera
	selectionne, il ne sera plus selectionnable, ainsi que les couples
	adjacents (8 au maximum).
	Les adresses correspondantes seront mises a 1.
  b.b- Calcul des minima (ou maxima):
       On boucle sur toutes les square distances du tableau TbDist2:
      - recherche d'un minimum (ou maximum) sur les echantillons,
      - calcul de l'extremum sur les courbes,
      - mise a jour du tableau TbSel.
-----------------------------------------------------------------------------*/
{
  myDone = Standard_False;

  const Handle(Cache)& aCache1 = myCache[0];
  const Handle(Cache)& aCache2 = myCache[1];
  Standard_NullObject_Raise_if ((aCache1.IsNull() || aCache2.IsNull()),
    "Extrema_GenExtCC::Perform()")

  Standard_Integer aNbU = aCache1->NbSamples(), aNbV = aCache2->NbSamples();
  Standard_OutOfRange_Raise_if ((aNbU < 2 ||aNbV < 2), "Extrema_GenExtCC::Perform()")

/*
a- Constitution du tableau des distances (TbDist2(0,NbU+1,0,NbV+1)):
   ---------------------------------------------------------------
*/

  //ensure that caches have been calculated
  if (!aCache1->IsValid())
    aCache1->CalculatePoints();
  if (!aCache2->IsValid())
    aCache2->CalculatePoints();

// Calcul des distances
  //initialization of variables in the same way as in CalculatePoints()
  Standard_Real PasU = aCache1->TrimLastParameter() - aCache1->TrimFirstParameter();
  Standard_Real PasV = aCache2->TrimLastParameter() - aCache2->TrimFirstParameter();
  Standard_Real U0 = PasU / aNbU / 100.;
  Standard_Real V0 = PasV / aNbV / 100.;
  PasU = (PasU - U0) / (aNbU - 1);
  PasV = (PasV - V0) / (aNbV - 1);
  U0 = aCache1->TrimFirstParameter() + (U0/2.);
  V0 = aCache2->TrimFirstParameter() + (V0/2.);

  const Handle(ArrayOfPnt)& aPntArray1 = aCache1->Points();
  const Handle(ArrayOfPnt)& aPntArray2 = aCache2->Points();
  Standard_Integer NoU, NoV;
  TColStd_Array2OfReal TbDist2(0, aNbU + 1, 0, aNbV + 1);
  for (NoU = 1; NoU <= aNbU; NoU++) {
    const Pnt& P1 = aPntArray1->Value (NoU);
    for (NoV = 1; NoV <= aNbV; NoV++) {
      const Pnt& P2 = aPntArray2->Value (NoV);
      TbDist2(NoU,NoV) = P1.SquareDistance(P2);
    }
  }

/*
b- Calcul des minima:
   -----------------
   b.a) Initialisations:
*/
//     - generales
  math_Vector Tol(1, 2);
  Tol(1) = myF.Tolerance();
  Tol(2) = myF.Tolerance();
  math_Vector UV(1,2), UVinf(1,2), UVsup(1,2);
  UVinf(1) = aCache1->TrimFirstParameter();
  UVinf(2) = aCache2->TrimFirstParameter();
  UVsup(1) = aCache1->TrimLastParameter();
  UVsup(2) = aCache2->TrimLastParameter();

//     - des 'bords' du tableau TbDist2
  for (NoV = 0; NoV <= aNbV+1; NoV++) {
    TbDist2(0,NoV) = RealLast();
    TbDist2(aNbU+1,NoV) = RealLast();
  }
  for (NoU = 1; NoU <= aNbU; NoU++) {
    TbDist2(NoU,0) = RealLast();
    TbDist2(NoU,aNbV+1) = RealLast();
  }

//     - du tableau TbSel(0,aNbU+1,0,aNbV+1) de selection des points
  TColStd_Array2OfInteger TbSel(0,aNbU+1,0,aNbV+1);
  TbSel.Init(0);

/*
   b.b) Calcul des minima:
*/
//     - recherche d un minimum sur la grille
  Standard_Integer Nu, Nv;
  Standard_Real Dist2;
  for (NoU = 1; NoU <= aNbU; NoU++) {
    for (NoV = 1; NoV <= aNbV; NoV++) {
      if (TbSel(NoU,NoV) == 0) {
	Dist2 = TbDist2(NoU,NoV);
	if ((TbDist2(NoU-1,NoV-1) >= Dist2) &&
	    (TbDist2(NoU-1,NoV  ) >= Dist2) &&
	    (TbDist2(NoU-1,NoV+1) >= Dist2) &&
	    (TbDist2(NoU  ,NoV-1) >= Dist2) &&
	    (TbDist2(NoU  ,NoV+1) >= Dist2) &&
	    (TbDist2(NoU+1,NoV-1) >= Dist2) &&
	    (TbDist2(NoU+1,NoV  ) >= Dist2) &&
	    (TbDist2(NoU+1,NoV+1) >= Dist2)) {

//     - calcul de l extremum sur la surface:
	  UV(1) = U0 + (NoU - 1) * PasU;
	  UV(2) = V0 + (NoV - 1) * PasV;


	  math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup);

      
//     - mise a jour du tableau TbSel 	  
	  for (Nu = NoU-1; Nu <= NoU+1; Nu++) {
	    for (Nv = NoV-1; Nv <= NoV+1; Nv++) {
	      TbSel(Nu,Nv) = 1;
	    }
	  }
	}
      } // if (TbSel(NoU,NoV)
    } // for (NoV = 1; ...
  } // for (NoU = 1; ...
/*
c- Calcul des maxima:
   -----------------
   c.a) Initialisations:
*/
//     - des 'bords' du tableau TbDist2
  for (NoV = 0; NoV <= aNbV+1; NoV++) {
    TbDist2(0,NoV) = RealFirst();
    TbDist2(aNbU+1,NoV) = RealFirst();
  }
  for (NoU = 1; NoU <= aNbU; NoU++) {
    TbDist2(NoU,0) = RealFirst();
    TbDist2(NoU,aNbV+1) = RealFirst();
  }

//     - du tableau TbSel(0,aNbU+1,0,aNbV+1) de selection des points
  TbSel.Init(0);
  /*for (NoU = 0; NoU <= aNbU+1; NoU++) {
    for (NoV = 0; NoV <= aNbV+1; NoV++) {
      TbSel(NoU,NoV) = 0;
    }
  }*/
/*
   c.b) Calcul des maxima:
*/
//     - recherche d un maximum sur la grille
  for (NoU = 1; NoU <= aNbU; NoU++) {
    for (NoV = 1; NoV <= aNbV; NoV++) {
      if (TbSel(NoU,NoV) == 0) {
	Dist2 = TbDist2(NoU,NoV);
	if ((TbDist2(NoU-1,NoV-1) <= Dist2) &&
	    (TbDist2(NoU-1,NoV  ) <= Dist2) &&
	    (TbDist2(NoU-1,NoV+1) <= Dist2) &&
	    (TbDist2(NoU  ,NoV-1) <= Dist2) &&
	    (TbDist2(NoU  ,NoV+1) <= Dist2) &&
	    (TbDist2(NoU+1,NoV-1) <= Dist2) &&
	    (TbDist2(NoU+1,NoV  ) <= Dist2) &&
	    (TbDist2(NoU+1,NoV+1) <= Dist2)) {

//     - calcul de l extremum sur la surface:
	  UV(1) = U0 + (NoU - 1) * PasU;
	  UV(2) = V0 + (NoV - 1) * PasV;

	  math_FunctionSetRoot S (myF,UV,Tol,UVinf,UVsup);
      
//     - mise a jour du tableau TbSel 	  
	  for (Nu = NoU-1; Nu <= NoU+1; Nu++) {
	    for (Nv = NoV-1; Nv <= NoV+1; Nv++) {
	      TbSel(Nu,Nv) = 1;
	    }
	  }
	}
      } // if (TbSel(NoU,NoV))
    } // for (NoV = 1; ...)
  } // for (NoU = 1; ...)
  myDone = Standard_True;
}
//=============================================================================

Standard_Boolean Extrema_GenExtCC::IsDone () const { return myDone; }
//=============================================================================

Standard_Integer Extrema_GenExtCC::NbExt () const
{
  StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()")
  return myF.NbExt();
}
//=============================================================================

Standard_Real Extrema_GenExtCC::SquareDistance (const Standard_Integer N) const
{
  StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
  Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
  return myF.SquareDistance(N);
}
//=============================================================================

void Extrema_GenExtCC::Points (const Standard_Integer N,
			    POnC& P1, POnC& P2) const
{
  StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
  Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
  myF.Points(N,P1,P2);
}