1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1313.c,v 1.12 2005-12-09 14:08:45 afr Exp $
45 *
46 */
47
48
49 #define S1313
50
51 #include "sislP.h"
52
53 /*
54 * Forward declarations.
55 * ---------------------
56 */
57 #if defined(SISLNEEDPROTOTYPES)
58 static void s1313_s9constline(SISLSurf *,double [],int,double,
59 SISLIntcurve *,int,int,int *);
60 #else
61 static void s1313_s9constline();
62 #endif
63
64 #if defined(SISLNEEDPROTOTYPES)
65 void
s1313(SISLSurf * ps1,double eimpli[],int ideg,double aepsco,double aepsge,double amax,SISLIntcurve * pintcr,int icur,int igraph,int * jstat)66 s1313(SISLSurf *ps1,double eimpli[],int ideg,double aepsco,double aepsge,
67 double amax,SISLIntcurve *pintcr,int icur,int igraph,int *jstat)
68 #else
69 void s1313(ps1,eimpli,ideg,aepsco,aepsge,amax,pintcr,icur,igraph,jstat)
70 SISLSurf *ps1;
71 double eimpli[];
72 int ideg;
73 double aepsco;
74 double aepsge;
75 double amax;
76 SISLIntcurve *pintcr;
77 int icur;
78 int igraph;
79 int *jstat;
80 #endif
81 /*
82 *********************************************************************
83 *
84 *********************************************************************
85 *
86 * PURPOSE : To march an intersection curve described by parameter pairs
87 * in the intersection curve, a B-spline surface and an
88 * implicit surface.
89 *
90 *
91 * INPUT : ps1 - Pointer to surface.
92 * eimpli - Description of the implicit surface
93 * ideg - The degree of the implicit surface
94 * ideg=1: Plane
95 * ideg=2; Quadric surface
96 * ideg=1001: Torus surface
97 * ideg=1003: Silhouette line parallel projection
98 * ideg=1004: Silhouette line perspective projection
99 * ideg=1005: Silhouette line circular projection
100 * aepsco - Not used.
101 * aepsge - Geometry resolution.
102 * amax - Not used.
103 * icur - Indicator telling if a 3-D curve is to be made
104 * 0 - Don't make 3-D curve
105 * 1 - Make 3-D curve
106 * 2 - Make 3-D curve and curves in parameter plane
107 * igraph - Indicator telling if the curve is to be outputted
108 * through function calls:
109 * 0 - don't output curve through function call
110 * 1 - output as straight line segments through
111 * s6move and s6line
112 *
113 *
114 *
115 * INPUT/OUTPUT:pintcr - The intersection curve. When comming as input
116 * only parameter values in the parameter plane
117 * exist. When comming as output the 3-D geometry
118 * and possibly the curve in the parameter plane
119 * of the surface are added.
120 * If the curves has already been generated in the
121 * topology part of the intersections, nothing will
122 * be done (i.e. not required). This will be the
123 * case when the intersection curve represents a
124 * constant parameter line in the parmeter plane
125 * of the surface.
126 *
127 * OUTPUT: jstat - status messages
128 * = 3 : Iteration stopped due to singular
129 * point or degenerate surface. A part
130 * of intersection curve may have been
131 * traced out. If no curve is traced out
132 * the curve pointers in the Intcurve
133 * object point to SISL_NULL.
134 * = 3 : Marching not succeded
135 * = 0 : ok
136 * < 0 : error
137 * = -185 : No points produced on intersection curve.
138 *
139 *
140 * METHOD :
141 * REFERENCES :
142 *
143 * HOW TO EXTEND THIS FUNCTION TO TREAT NEW PROBLEMS.
144 *
145 * This function is built as a general function for treating the combination
146 * of an implicit description and a parametric surface, when introducing
147 * new problems to be solved this structure can be utilized:
148 *
149 * 1. Define the implicit degree of the problem. If it is a special problem
150 * utilize ideg>1000.
151 * 2. Determine how many derivatives have to be calculated to support the
152 * function. For ideg=1,2,1001 derivatives up two have been calculated.
153 * for ideg=1003,1004,1005, derivatives up to and including 3 have been calculated.
154 * The number of derivatives influence the local variables kder, ksize,
155 * ksizem3 that are found in the functions:
156 * s1306,s1309,s1313,s1331,s9clipimp,s9iterimp,s9adsimp,s9boundimp
157 *
158 * 3. The function s1309 branch on itype to decide the distance in the
159 * current iteration step. This function has to be updated to support
160 * the new type of problem. For ideg=1,2 and 1001 this function
161 * currently calculates distance. For ideg=1003,1004,1005 it calculates an angle.
162 *
163 * 4. The function s1331 that calculates derivatives of the combination of
164 * the implicit problem and the surface have to support the new problem
165 *
166 * 5. Look also into s9iterimp and s9boundimp to update the convergence
167 * branching on problem type.
168 *
169 *
170 *-
171 * CALLS :
172 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway, 2. July 1988
173 * Revised by : Tor Dokken, SI, Oslo, Norway, 22. January 1988,
174 * Test for degeneracy and singularities included
175 * Revised by : Tor Dokken, SI, Oslo, Norway, March 1989
176 * Automatic generation of maximal step length, improved
177 * marching close to singularities
178 * Revised by : Correction of error testing
179 * Revised by : Mike Floater, SI, 1991-01
180 * Improved the routine for parallel silhouettes (ideg=1003) and
181 * added perspective and circular silhouettes (ideg=1004,ideg=1005)
182 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, Dec. 1994. Added check for
183 * SISL_NULL 'pintcr' and to avoid re-generating the geometry when it has
184 * already been generated in the topology part (constant curve, type 9).
185 * This fixes memory problems.
186 *
187 *********************************************************************
188 */
189 {
190 int ki,kj,kl; /* Control variables in for loops */
191 int kcont; /* Stop condition for loop */
192 int kk,kn; /* Dummy variables */
193 int kstpch; /* Status of iteration step */
194 int kpoint; /* Number of points in guide curve */
195 int kpar1; /* Number of parameter direction in 1st. obj */
196 int kpar2; /* Number of parameter direction in 2st. obj */
197 int kpar; /* Indicater tellin if S1359 shall make
198 parametrization or use parametrization
199 in spar */
200 int ktype; /* Type of intersection curve */
201 int klfu=0; /* Pointers into knot vectors */
202 int klfv=0; /* Pointers into knot vectors */
203 int kder = 2; /* Calculate up to second derivatives */
204 int kdim = 3; /* The dimension of the space we work in */
205 int kfirst = 0; /* Indicator telling if first guide point degenerate */
206 int klast = 0; /* Indicator telling if last guide point degenerate */
207 int kpos = 0; /* Position of error */
208 int kstat,kstat1; /* Status variable returned form routine */
209 int kmaxinf=0; /* Number of entries object that can be stored
210 in s3dinf, sp1inf, sp2inf */
211 int knbinf=0; /* Number of entries stored so far on s3dinf,
212 sp1inf and sp2inf */
213 int kstart; /* Start point for iteration among guide pnts*/
214 int kguide; /* Current guide point */
215 int kdir; /* March direction */
216 int kgdir; /* Direction we march guide point vector */
217 int krem,krem1,krem2; /* Remember if we step in or out of patch */
218 int kbound; /* Dummy variiable */
219 int koutside_resolution; /* Flag telling if current seg. outside res. */
220 int ksize; /* Number of doubles for storage of derivateves
221 and normal vector */
222 int ksizem3; /* ksize - 3 */
223 int kdiv=0; /* Flag remembering if iteration diverged */
224 int knb1=0; /* Remember number of points after marching
225 in first marching direction */
226 int kgd1=0; /* Remeber last guide point used in first
227 marching direction */
228 double *scorpnt=SISL_NULL; /* Corrected marching points */
229 double *scorpar=SISL_NULL; /* Corrected marching parameter values */
230 double smidd[6]; /* Description of midpoint and tangent of
231 current Bezier segment */
232 double tcurstep; /* Current step length */
233 double tdist; /* Error at middle of current Bezier segement*/
234 double tang; /* Angle error at midpoint Bezier segement */
235 double tnew; /* Candidate for new step length */
236 double tfak; /* How much is the step length to be reduced */
237 double *start; /* Pointer to start of current segment */
238 double *st; /* Pointer to knot vector */
239 double sval1[2]; /* Limits of parameter plane in first SISLdir */
240 double tref1,tref2; /* Reference values for knot vectors */
241 double sval2[2]; /* Limits of parameter plane in second SISLdir */
242 double tstep; /* Iteration step length */
243 double tmax; /* Local maximal step length */
244 double tstartstp; /* Start step length */
245 double trad; /* Radius of curvature */
246 double tval[6]; /* Dummy array in s1331 */
247 double *spar=SISL_NULL; /* Parametrization of points in Hermit interp*/
248 double spar1[2]; /* Parameter pair of current point surface 1 */
249 double spar2[2]; /* Parameter pair of boundarypoint surface 1 */
250 double siparmid[2]; /* Parameter value at middle of Bezier segment*/
251 double sipar1[2]; /* Parameter pair iteration point surface 1 */
252 double simiddpnt[10]; /* Middle point and tangent of segment */
253 double simiddpar[7]; /* Parameter value at middle point of segment*/
254 double *sgpar1=SISL_NULL; /* Parameter pairs of guide point in surf 1 */
255 double *sgpar2=SISL_NULL; /* Parameter pairs of guide point in surf 2 */
256 double *sgpar=SISL_NULL; /* guide points used */
257 double *sgd1 = SISL_NULL; /* 0-2 derivative of guide point + normal
258 of first object */
259 double spnt1[33]; /* Info on current point in first surface */
260 double spnt2[33]; /* Info on boundary point in first surface */
261 double sipnt1[33]; /* Info on iteration point in first surface */
262 /* For spnt1, sipnt1, the */
263 /* information is stored 3-tuppels in the */
264 /* following sequence */
265 /* Position, (1,0)-der, (0,1)-der, */
266 /* (2,0)-der, (1,1)-der, (0,2)-der and normal*/
267 /* This is compatible with output of s1421 */
268 double snorm1[3]; /* Normal vector of implicit surface */
269 double snorm2[3]; /* Normal vector of implicit surface */
270 double startg[3]; /* Start tangent of iteration */
271
272 double *snxt1; /* SISLPoint in ps1 we have accepted */
273 double *snxp1; /* Parameter value belonging to snxt1 */
274 double *s3dinf=SISL_NULL; /* Pointer to array used for storing 3-D position
275 tangent, curvature and radius of curvature found
276 during the marching process if possible */
277 double *sp1inf=SISL_NULL; /* Pointer to array used for storing position
278 tangent, curvature and radius of curvature found
279 in the first parameter plane during the
280 marching process */
281 double start1[33]; /* Description of start point in ps1 */
282 double stpar1[2]; /* Parameter pair belonging to start1 */
283 double sdum1[3],sdum2[3];/* Dummy vectors */
284 double tdum,tdump; /* Dummy variable */
285 double *sp1=SISL_NULL; /* Pointer used when moving information */
286 double *sp2=SISL_NULL; /* Pointer used when moving information */
287 double stdum[10]; /* Dummy array used when moving information */
288 double *stang; /* Pointer to tangent of current point */
289 double *sptang; /* Pointer to tangent in parameter plane */
290 double *spoint; /* Pointer to current point */
291 double t1distgd,t2distgd;/* Distances to guide points */
292 double tlnorm; /* Length of normal vector */
293 double tltan1,tltan2; /* Tangent lengths */
294 SISLCurve *q3dcur=SISL_NULL;/* Pointer to 3-D curve */
295 SISLCurve *qp1cur=SISL_NULL;/* Pointer to curve in first parameter plane*/
296 double sdiffcur[3]; /* Difference between current and previous point found */
297 double sdiffprev[3]; /* Difference between previous point and the one before that */
298
299
300 *jstat = 0;
301
302 if (pintcr == SISL_NULL) goto err150;
303
304
305 /* Check if the geometry already has been generated in the topology part.
306 This will be the case if the geometry is along a constant parameter line. */
307
308 if (pintcr->itype == 9) goto out;
309
310
311
312 /* Make maximal step length based on box-size of surface */
313
314 sh1992su(ps1,0,aepsge,&kstat);
315 if (kstat < 0) goto error;
316
317 tmax = MAX(ps1->pbox->e2max[0][0] - ps1->pbox->e2min[0][0],
318 ps1->pbox->e2max[0][1] - ps1->pbox->e2min[0][1]);
319 tmax = MAX(tmax,ps1->pbox->e2max[0][2] - ps1->pbox->e2min[0][2]);
320
321 if (amax>DZERO) tmax = MIN(tmax,amax);
322
323
324 /* If ideg=1,2 or 1001 then only derivatives up to second order
325 are calculated, then 18 doubles for derivatives and 3 for the
326 normal vector are to be used for calculation of points in the
327 spline surface. For ideg=1003,1004,1005 we have a silhouette curve and
328 derivatives up to the third are to be calculated,
329 thus 30 +3 a total of 33 doubles are to be calculated */
330
331 if (ideg==1003 || ideg==1004 || ideg==1005)
332 {
333 kder = 3;
334 ksize = 33;
335 }
336 else
337 {
338 ksize = 21;
339 kder =2;
340 }
341 ksizem3 = ksize -3;
342
343
344 /* Find a none singular start point for the marching process */
345
346 kpoint = pintcr->ipoint;
347 kpar1 = pintcr->ipar1;
348 kpar2 = pintcr->ipar2;
349 sgpar1 = pintcr->epar1;
350 sgpar2 = pintcr->epar2;
351 ktype = pintcr->itype;
352
353
354 /* Initiate pointers to intersection curve and intersection curve in
355 parameter plane */
356
357 pintcr -> pgeom = SISL_NULL;
358 pintcr -> ppar1 = SISL_NULL;
359 pintcr -> ppar2 = SISL_NULL;
360
361
362 /* Initiate parameter direction boundaries */
363 kk = ps1 -> ik1;
364 kn = ps1 -> in1;
365 st = ps1 -> et1;
366 sval1[0] = st[kk-1];
367 sval1[1] = st[kn];
368 tref1 = (double)3.0*MAX(fabs(*sval1),fabs(*(sval1+1)));
369 kk = ps1 -> ik2;
370 kn = ps1 -> in2;
371 st = ps1 -> et2;
372 sval2[0] = st[kk-1];
373 sval2[1] = st[kn];
374 tref2 = (double)3.0*MAX(fabs(*sval2),fabs(*(sval2+1)));
375
376 /* Test that first object has 2 parameter direction and second object 0*/
377
378 if (kpar1 == 2 && kpar2 == 0)
379 {
380 /* Everithing is ok */
381 ;
382 }
383 else if (kpar1 == 0 && kpar2 == 2)
384 {
385 sgpar1 = sgpar2;
386 }
387 else
388 {
389 goto err123;
390 }
391
392 /* To support closed curve the first guide point must be copied after
393 the last guide point */
394
395 if((sgpar=newarray(2*kpoint+2,DOUBLE)) == SISL_NULL) goto err101;
396 memcopy(sgpar,sgpar1,2*kpoint,DOUBLE);
397 if (ktype ==2 || ktype == 3)
398 {
399 /* Closed curve copy first guide point to end of string of guide points */
400 memcopy(sgpar+2*kpoint,sgpar1,2,DOUBLE);
401 kpoint = kpoint + 1;
402 }
403
404 /* THE POINTS , TANGENT, CURVATURE AND RADIUS OF CURVATURE FOUND DURING
405 * THE MARCHING PROCESS SHOULD ALL BE STORED IN ARRAYS. ALLOCATE ONE ARRAY
406 * FOR 3-D INFORMATION AND ONE ARRAY FOR INFORMATION IN THE PARAMETER PLANE.
407 * THESE ARRAYS ARE GIVEN AN INITIAL CAPACITY OF STORING 100 POINTS WITH
408 * OTHER INFORMATION.
409 * IF THEY ARE TOO SHORT THEY WILL BE REALLOCATED AT A LATER STAGE.
410 *
411 * SINCE THE MARCHING WILL GO IN BOTH DIRECTIONS WE WILL HAVE TO TURN THE
412 * INFORMATION FOUND WHEN MARCHING IN NEGATIVE DIRECTION, SO THAT IT CAN
413 * BE COMBINED WITH THE INFORMATION FOUND WHEN WE ARE MARCHING IN POSITVE
414 * DIRECTION.
415 */
416
417 kmaxinf = 100;
418 s3dinf = newarray(10*kmaxinf,DOUBLE);
419 if (s3dinf == SISL_NULL) goto err101;
420 sp1inf = newarray(7*kmaxinf,DOUBLE);
421 if (sp1inf == SISL_NULL) goto err101;
422
423 /* Evaluate 0-1-2nd. derivative + normal of all guide points in the surface,
424 first allocate arrays for storing the information, check that the points
425 have a defined normal, and that the combination of the implicit surface
426 and the surface defines a tangent direction in the curve */
427
428 sgd1 = newarray(ksize*kpoint,DOUBLE);
429 if (sgd1==SISL_NULL) goto err101;
430
431 kpos = 5;
432
433 /* Initiate kstart to point at no point */
434
435 kstart = 0;
436
437 for (ki=0,kj=0,kl=0 ; ki<kpoint ; ki++,kj+=2,kl+=ksize)
438 {
439 s1421(ps1,kder,&sgpar[kj],&klfu,&klfv,&sgd1[kl],&sgd1[kl+ksizem3],&kstat);
440 if (kstat<0) goto error;
441
442
443
444 if (ideg == 1003 || ideg == 1004 || ideg == 1005)
445 {
446 /* Find length of normal vector and tangent vectors */
447
448 tlnorm = s6length(&sgd1[kl+ksizem3],kdim,&kstat);
449 tltan1 = s6length(&sgd1[kl+kdim],kdim,&kstat);
450 tltan2 = s6length(&sgd1[kl+kdim+kdim],kdim,&kstat);
451
452 /* The cross product satisifes the follwing conditions:
453 length(axb) = length(a) length(b) sin(angle(a,b)).
454 Thus the angle between the two vectors can be found, close to 0
455 sin(a) is a good approximation of a
456 */
457 if (tlnorm == (double)0.0 || tltan1 ==(double)0.0 || tltan2 == (double)0.0)
458 tang = (double)0.0;
459 else
460 tang = tlnorm/(tltan1*tltan2);
461
462 /* Silhouette line calculation no normal can be found in implicit surface
463 accept point as candidate start point if tang greater tha angular
464 resolution */
465
466 if (tang >= ANGULAR_TOLERANCE) kstart = ki+1;
467 }
468 else
469 {
470 /* Make direction of intersection curve */
471
472 s1306(&sgd1[kl],&sgpar[kj],eimpli,ideg,s3dinf,sp1inf,&kstat);
473 if (kstat < 0) goto error;
474
475
476 /* Remember if start, internal or end point */
477
478 if (kstat != 2)
479 {
480 if (ki == 0) kfirst = 1;
481 else if (ki == kpoint-1) klast = kpoint;
482 else kstart = ki+1;
483 }
484 }
485 }
486
487 /* Check if only degenerate points or singularities exist on the
488 intersection curve */
489
490 if (kstart == 0)
491 {
492 /* No internal nondegenerate point exits, start marching from first
493 or last point if possible */
494 if (kfirst != 0 && ktype != 5 && ktype != 7) kstart = kfirst;
495 else if (klast != 0 && ktype != 6 &&
496 ktype != 7 && ktype != 3) kstart = klast;
497 else if (kfirst != 0) kstart = kfirst;
498 else if (klast != 0) kstart = klast;
499 else goto interpolate;
500 }
501
502 /* To speed up the marching process when many guide points are given,
503 remove guide points that are not at the start, end or the start point */
504
505 if (kpoint >2 && (kstart==1 || kstart==kpoint) )
506 {
507 /* No internal guide point necessary, copy last point
508 to second point */
509
510 memcopy(sgd1+ksize,sgd1+ksize*(kpoint-1),ksize,DOUBLE);
511 memcopy(sgpar+2,sgpar+2*(kpoint-1),2,DOUBLE);
512
513 if (kstart == kpoint) kstart = 2;
514 kpoint = 2;
515 }
516 else if (kpoint>2)
517 {
518 /*Internal guide point exists, copy this to second position and
519 copy end point to third position */
520
521 memcopy(sgd1+ksize,sgd1+ksize*(kstart-1),ksize,DOUBLE);
522 memcopy(sgpar+2,sgpar+2*(kstart-1),2,DOUBLE);
523
524 memcopy(sgd1+2*ksize,sgd1+ksize*(kpoint-1),ksize,DOUBLE);
525 memcopy(sgpar+4,sgpar+2*(kpoint-1),2,DOUBLE);
526
527 kpoint = 3;
528 kstart = 2;
529 }
530
531
532
533 /* Remember description of start point in both surfaces,
534 * copy point indicated by kstart into spnt1,spar1 */
535
536 memcopy(spnt1,sgd1+ksize*(kstart-1),ksize,DOUBLE);
537 memcopy(spar1,sgpar+2*(kstart-1),2,DOUBLE);
538
539 /* Make position, unit tangent, curvature and radius of curvature for
540 * start point of iteration, store them in the arrays just allocated */
541
542 kpos = 10;
543 s1306(spnt1,spar1,eimpli,ideg,s3dinf,sp1inf,&kstat);
544
545 if (kstat<0) goto error;
546
547 /* Test if singular point reached */
548
549 if (kstat == 2) goto war03;
550
551 /* Remember start tangent */
552 memcopy(startg,s3dinf+3,3,DOUBLE);
553
554
555 /* Iterate intersection point down to the intersection curve */
556
557 tstep = DZERO;
558 s9iterimp(s3dinf,spnt1,spar1,ps1,eimpli,ideg,tstep,
559 aepsge,sipnt1,sipar1,&kstat);
560 if (kstat < 0) goto error;
561
562 if (kstat==0 && s6dist(spnt1,sipnt1,3) > aepsge)
563 {
564
565 /* Nonsingular point found and adjustment greater than aepsge,
566 copy result of iteration into spnt1,spar1 */
567
568 memcopy(spnt1,sipnt1,ksize,DOUBLE);
569 memcopy(spar1,sipar1,2,DOUBLE);
570 }
571
572 /* Remember start point */
573
574 if (kstat==0)
575 {
576 memcopy(start1,sipnt1,ksize,DOUBLE);
577 memcopy(stpar1,sipar1,2,DOUBLE);
578 }
579 else
580 {
581 memcopy(start1,spnt1,ksize,DOUBLE);
582 memcopy(stpar1,spar1,2,DOUBLE);
583 }
584
585 /* Make position, unit tangent, curvature and radius of curvature for
586 start point of iteration, store them in the arrays just allocated */
587
588 kpos = 10;
589 s1306(start1,stpar1,eimpli,ideg,s3dinf,sp1inf,&kstat);
590
591 if (kstat<0) goto error;
592
593
594 /* Test if singular point reached */
595
596 if (kstat == 2) goto war03;
597
598 /* Remember that start point is already stored */
599
600 knbinf = 1;
601
602 /* Make step length based on 3-D radius of curvature, tolerances and
603 maks step length */
604
605 kpos = 20;
606 tstep = s1311(s3dinf[9],aepsge,tmax,&kstat);
607 if (kstat<0) goto error;
608 tstartstp = tstep;
609
610 /* STEP IN BOTH DIRECTIONS FROM THE FOUND START POINT */
611
612 /* Indicate that direction in guide point array not determined */
613 kguide = kstart;
614 kgdir = 0;
615
616 for (kdir=1;kdir<3;kdir++)
617 {
618
619 if (kdir == 2)
620 {
621
622 /* Remember result of marching in first direction */
623
624 knb1 = knbinf;
625 kgd1 = kguide;
626
627
628 /* If the previous step direction made no points then knbinf==0. To
629 enable the marching we start from the same start point as the
630 previous step direction, thus in this case knbinf should be 1. */
631
632 knbinf = MAX(1,knbinf);
633
634 /* We now step in the second step direction. Turn the sequence of
635 the points found as well as change tangent directions */
636
637 /* First interchange 3-D info */
638
639 for (sp1=s3dinf,sp2=s3dinf+10*(knbinf-1) ; sp1<sp2 ; sp1+=10,sp2-=10)
640 {
641 memcopy(stdum,sp1, 10,DOUBLE);
642 memcopy(sp1 ,sp2, 10,DOUBLE);
643 memcopy(sp2 ,stdum,10,DOUBLE);
644 }
645
646 for (sp1=s3dinf+3;sp1<s3dinf+10*knbinf;sp1+=10)
647 {
648 sp1[0] = - sp1[0];
649 sp1[1] = - sp1[1];
650 sp1[2] = - sp1[2];
651 }
652
653 /* Then interchange info in parameter plane */
654
655 for (sp1=sp1inf,sp2=sp1inf+7*(knbinf-1) ; sp1<sp2 ; sp1+=7,sp2-=7)
656 {
657 memcopy(stdum,sp1 ,7,DOUBLE);
658 memcopy(sp1 ,sp2 ,7,DOUBLE);
659 memcopy(sp2 ,stdum,7,DOUBLE);
660 }
661
662 for (sp1=sp1inf+2;sp1<sp1inf+7*knbinf;sp1+=7)
663 {
664 sp1[0] = - sp1[0];
665 sp1[1] = - sp1[1];
666 /* sp1[2] = - sp1[2]; */
667 }
668
669 /* Turn direction of remembered start tangent */
670
671 for (ki=0;ki<3;ki++)
672 startg[ki] = -startg[ki];
673
674
675 /* Update spnt1 and spar1 to have the start point values */
676
677 memcopy(spnt1,start1,ksize,DOUBLE);
678 memcopy(spar1,stpar1,2,DOUBLE);
679
680 /* Turn the direction we march the guide point vector, and set current
681 guide point to kstart */
682
683 kgdir = -kgdir;
684 kguide = kstart;
685
686 /* Update step length */
687 tstep = tstartstp;
688 }
689
690
691 stang = s3dinf + 10*(knbinf-1) + 3;
692 kpos = 30;
693
694 /* Step direction ok, perform marching until stop condition reached */
695
696 kcont = 1;
697
698 while (kcont)
699 {
700
701 /* We must make sure that we are not stepping past a guide point.
702 * Thus if we get close to a guide point, make sure that we step
703 * through this. The direction we travers the guide point array
704 * might not have been determined yet. Thus we have to test in
705 * both directions in guide point array.
706 *
707 * Remember how we step in the varaible kstpch:
708 * kstpch = -1 : Try to step to previous guide point
709 * kstpch = 0 : Try not to step through guide point
710 * kstpch = 1 : Try to step to next guide point
711 * kstpch = 3 : Step to start point and stop marching
712 * kstpch = 4 : Don't step through guide point, candidate
713 * end point of segement found in iteration loop
714 */
715
716
717 kstpch = 0;
718 stang = s3dinf + 10*(knbinf-1) + 3;
719 sptang = sp1inf + 7*(knbinf-1) + 2;
720 if (kgdir >=0)
721 {
722
723 /* We are stepping in positive direction in guide point vector
724 * calculate distance to next guide point. If the guide point
725 * is lying closer than the step length to the current point
726 * we should step directly to this point provided that the cross
727 * product of the normal vectors at current point and at the
728 * guide point have the same direction, e.g. that their scalar
729 * product is positiv */
730
731 t1distgd = (double)2.0*tstep;
732 if (kguide < kpoint)
733 {
734 /* Decide if we should step through the guide point */
735
736 kpos = 40;
737 t1distgd = s9adsimp(spnt1,spar1,eimpli,ideg,
738 &sgd1[kguide*ksize],
739 &sgpar[kguide*2],
740 stang,sptang,tstep,&kstat);
741 if (kstat < 0) goto error;
742 if (kstat == 1)
743 {
744 /* Step through guide point remember this */
745
746 kstpch = 1;
747 snxt1 = sgd1 + ksize*kguide;
748 snxp1 = sgpar + 2*kguide;
749 tstep = MIN(tstep,t1distgd);
750 }
751 }
752 }
753
754 if (kgdir <=0)
755 {
756
757 /* We are stepping in negative direction in guide point vector
758 * calculate distance to previous guide point. If the guide point
759 * is lying closer than the step length to the current point
760 * we should step directly to this point provided that the cross
761 * product of the normal vectors at current point and at the
762 * guide point have the same direction, e.g. that their scalar
763 * product is positiv */
764
765 if (1 < kguide)
766 {
767 /* Decide if we should step through the guide point */
768
769 kpos = 50;
770 t2distgd = s9adsimp(spnt1,spar1,eimpli,ideg,
771 &sgd1[(kguide-2)*ksize],&sgpar[2*kguide-4],
772 stang,sptang,tstep,&kstat);
773 if (kstat<0) goto error;
774 if ((kstat == 1 &&kstpch == 0) ||
775 (kstat == 1 && kstpch == 1 && t2distgd < t1distgd))
776 {
777 /* Step through guide point remember this */
778
779 kstpch = -1;
780 snxt1 = sgd1 + ksize*(kguide-2);
781 snxp1 = sgpar + 2*(kguide-2);
782 tstep = MIN(tstep,t2distgd);
783 }
784 }
785 }
786
787 /* Check if we step through the start point. Should only be done if at least
788 three points found in this marching direction, a full closed curve will
789 require 6 segments. */
790 if ((kdir==1 && knbinf>3) || (kdir==2 && knbinf>knb1+3))
791 {
792
793 kpos = 60;
794 tdum = s9adsimp(spnt1,spar1,eimpli,ideg,start1,stpar1,stang,sptang,
795 tstep,&kstat);
796 if (kstat < 0) goto error;
797 if (kstat == 1)
798 {
799 /* Step to start point remember this */
800
801 kstpch = 3;
802 snxt1 = start1;
803 snxp1 = stpar1;
804 tstep = MIN(tstep,tdum);
805 }
806 }
807
808
809 /* At this stage kstpch=0 if we have not reached a guide point or
810 if we have not reached the start point of the iteration.
811
812 Now we want to find a Bezier segement that is lying within the
813 geometric tolerance that is approximating the intersection curve.
814 If a guide point is reached (kstpch=-1 or 1), then we have a
815 candidate for the end point of the Bezier segement. If the start
816 point is reached (kstpch=3) then we also have a candidate end point
817 for the segment.
818
819 The next loop use kstpch to indicate if we have a candidate
820 end point for the segment:
821
822 kstpch==0 : No candidate end point exists
823 kstpch!=0 : Candidate end point exists
824
825
826 To indicate if the segement is within the resolution we
827 use koutside_resolution:
828
829 koutside_resolution==0 : Segment outside resolution
830 koutside_resolution!=0 : Segment inside resolution
831 */
832
833 koutside_resolution = 0;
834
835 /* Make sure that there is enough space for one more point */
836 if (knbinf>=kmaxinf)
837 {
838 kmaxinf = kmaxinf + 100;
839 s3dinf = increasearray(s3dinf,((3*kdim+1)*kmaxinf),DOUBLE);
840 if (s3dinf==SISL_NULL) goto err101;
841 sp1inf = increasearray(sp1inf,7*kmaxinf,DOUBLE);
842 if (sp1inf==SISL_NULL) goto err101;
843 }
844
845 /* Make description of candidate endpoint if it exists and store it */
846
847 if (kstpch != 0)
848 {
849 s1306(snxt1,snxp1,eimpli,ideg,s3dinf+10*knbinf,
850 sp1inf+7*knbinf,&kstat);
851 if (kstat<0) goto error;
852
853 /* It is allowed to jump on to a singular point */
854
855 /* Make sure that the tangents of previous and the
856 candidate end point point in the same direction */
857
858 if (knbinf>0)
859 tdum = s6scpr(s3dinf+10*(knbinf-1)+3,
860 s3dinf+10*knbinf+3,kdim);
861 else
862 tdum = s6scpr(startg,s3dinf+3,kdim);
863
864
865 if (tdum < DZERO)
866 {
867 /* Change tangent direction 3-D and in parameter plane */
868 sp1 = s3dinf + 10*knbinf + 3;
869 sp1[0] = -sp1[0];
870 sp1[1] = -sp1[1];
871 sp1[2] = -sp1[2];
872 sp1 = sp1inf + 7*knbinf + 2;
873 sp1[0] = -sp1[0];
874 sp1[1] = -sp1[1];
875 }
876
877 /* Copy the candidate point to spnt2 and spar2 */
878
879 memcopy(spnt2,snxt1,ksize,DOUBLE);
880 memcopy(spar2,snxp1,2,DOUBLE);
881 }
882
883 while (kstpch == 0 || koutside_resolution == 0)
884 {
885 spoint = s3dinf + 10*(knbinf-1); // To avoid confusing valgrind
886 if (kstpch!=0)
887 {
888 /* Candidate end point exist, iterate to find point close
889 to the midpoint of the Bezier segement */
890
891
892 /* Decide if Hermit shape acceptable and find position and
893 tangent at midpoint of segment */
894
895 start = s3dinf + 10*(knbinf-1);
896
897 s1361(start,start+10,3,smidd,smidd+3,&kstat);
898 if (kstat<0) goto error;
899
900 tcurstep = DZERO;
901 spoint = smidd;
902 }
903 else
904 {
905
906 /* Iterate to find end point of segment */
907
908 /* ITERATE by intersecting the two surface and the plane
909 * defined by current point (s3dinf), the tangent (s3dinf+3)
910 * and the step length */
911
912 spoint = s3dinf + 10*(knbinf-1);
913 tcurstep = tstep;
914 }
915
916 /* Perform the actual iteration */
917
918 kpos = 70;
919 s9iterimp(spoint,spnt1,spar1,ps1,eimpli,ideg,tcurstep,
920 aepsge,sipnt1,sipar1,&kstat);
921 if (kstat < 0) goto error;
922
923 /* Initiate distance between midpoint and iteration point
924 to -1 to enable detection of divergence */
925
926 tdist = -1;
927
928 /* Check if iteration has converged */
929
930 if (kstat == 2)
931 {
932 /* Iteration has diverged, half step length if possible,
933 find new endpoint of segement. */
934
935 kstpch = 0;
936 koutside_resolution = 0;
937 }
938 else if(kstat == 1 && kstpch != 0)
939 {
940 /* The point found is closer to the input point than
941 the relative computer resolution or is a singular point.
942 Half step length if possible, find new endpoint of
943 segement. */
944
945 kstpch = 0;
946 koutside_resolution = 0;
947 }
948 else if (kstpch!=0)
949 {
950
951
952 /* Make description of intersection point */
953
954 s1306(sipnt1,sipar1,eimpli,ideg,simiddpnt,simiddpar,&kstat);
955 if (kstat<0) goto error;
956
957 if (kstat != 2)
958 {
959 /* We iterated to find midpoint of segment, test if
960 it is within resolution */
961
962 tdist = s6dist(simiddpnt,smidd,3);
963 tang = s6ang(simiddpnt+3,smidd+3,3);
964 }
965
966 /* If point is singular or not within resolution a new
967 Hermit segment has to be made */
968
969 if (kstat == 2 || (fabs(tdist) > aepsge ||
970 (fabs(tang) > ANGULAR_TOLERANCE && tstep>aepsge)))
971 {
972 kstpch = 0;
973 koutside_resolution = 0;
974 }
975 else
976 {
977 /* Segment within tolerance.
978 * Check that the relationship between the two surfaces
979 * has not been interchanged, by making the cross product
980 * of the normal vectors in current point and the point
981 * found by iteration. Then make the scalar product of
982 * these vectors. If the scalar product is negative then
983 * we have either jumped to another branch or passed a
984 * singularity, iterprete this as the iteration has diverged
985 * In addition we don't want the direction of the tangents
986 * change to much. We set a limit of approximately 41 degrees
987 * by testing on a cosin value of 0.75
988 * Make normal vectors in implicit surface for both points
989 * Make also sure that the curve in the parameter plane
990 * does not turn more than 90 degrees.
991 */
992
993 s1331(spnt1,eimpli,ideg,-1,tval,snorm1,&kstat);
994 if (kstat < 0) goto error;
995 s1331(spnt2,eimpli,ideg,-1,tval,snorm2,&kstat);
996 if (kstat < 0) goto error;
997
998 if(ideg == 1003 || ideg == 1004 || ideg == 1005)
999 {
1000 tdum = s6scpr(snorm1,snorm2,kdim);
1001 }
1002 else
1003 {
1004 s6crss(spnt1+ksizem3,snorm1,sdum1);
1005 (void)s6norm(sdum1,kdim,sdum1,&kstat);
1006 if (kstat < 0) goto error;
1007 s6crss(sipnt1+ksizem3,snorm2,sdum2);
1008 (void)s6norm(sdum2,kdim,sdum2,&kstat);
1009 if (kstat < 0) goto error;
1010 tdum = s6scpr(sdum1,sdum2,kdim);
1011 }
1012
1013 s6diff(spar2,spar1,2,sdum1);
1014 tdump = s6scpr(sdum1,sp1inf+7*(knbinf-1)+2,2);
1015 /* The test below was added to detect the case when the
1016 normal of the parametric surface turns */
1017 if(s6scpr(spnt1+ksizem3,sipnt1+ksizem3,kdim) < 0.0)
1018 {
1019 tdum = -tdum;
1020 }
1021 if (tdum == DZERO)
1022 {
1023 double tl1,tl2;
1024
1025 /* If one of the tangents have zero length, accept
1026 segment */
1027
1028 tl1 = s6length(sdum1,kdim,&kstat);
1029 tl2 = s6length(sdum2,kdim,&kstat);
1030
1031 if (tl1 == DZERO || tl2 == DZERO)
1032 koutside_resolution = 1;
1033 else
1034 {
1035 /* Find new end point of segment */
1036 koutside_resolution = 0;
1037 kstpch = 0;
1038 }
1039
1040 }
1041 /*else if (tdum <= (double)0.75 || tdump <= (double)0.0)*/
1042 else if (tdum <= (double)0.75 || tdump <= -REL_PAR_RES)
1043 {
1044 /* Find new end point of segment */
1045 koutside_resolution = 0;
1046 kstpch = 0;
1047 }
1048 else
1049 {
1050 koutside_resolution = 1;
1051 }
1052 }
1053 }
1054 else
1055 {
1056 /* We iterated to find end point of segment, update pointer */
1057
1058 memcopy(spnt2,sipnt1,ksize,DOUBLE);
1059 memcopy(spar2,sipar1,2,DOUBLE);
1060
1061 s1306(sipnt1,sipar1,eimpli,ideg,s3dinf+10*knbinf,
1062 sp1inf+7*knbinf,&kstat);
1063 if (kstat<0) goto error;
1064
1065 /* Make sure that the tangents of previous and the new point
1066 point in the same direction, singular end point allowed */
1067
1068 if (knbinf>0)
1069 tdum = s6scpr(s3dinf+10*(knbinf-1)+3,
1070 s3dinf+10*knbinf+3,kdim);
1071 else
1072 tdum = s6scpr(startg,s3dinf+3,kdim);
1073
1074
1075 if (tdum < DZERO)
1076 {
1077 /* Change tangent direction 3-D and in parameter plane */
1078 sp1 = s3dinf + 10*knbinf + 3;
1079 sp1[0] = -sp1[0];
1080 sp1[1] = -sp1[1];
1081 sp1[2] = -sp1[2];
1082 sp1 = sp1inf + 7*knbinf + 2;
1083 sp1[0] = -sp1[0];
1084 sp1[1] = -sp1[1];
1085 }
1086
1087 /* Indicate that end point accepted */
1088
1089 kstpch = 4;
1090 koutside_resolution = 0;
1091 }
1092
1093 /* If the segment is accepted. check if we cross the
1094 boundary of the patch */
1095
1096 if (kstpch !=0 && koutside_resolution == 1)
1097 {
1098
1099 /* Check if curve between start and endpoint cross the
1100 boundary */
1101
1102 memcopy(siparmid,sipar1,2,double);
1103
1104 s1305(spar1,spar2,sval1,sval2,&kbound,sipar1,&kstat);
1105 if (kstat<0) goto error;
1106
1107
1108 if(kstat==0 || kstat==4)
1109 {
1110 /* Set pointer to tangents at start point */
1111 ki = 7*(knbinf-1)+2;
1112
1113 if( (DEQUAL(spar1[1]+tref2,sval2[0]+tref2) && sp1inf[ki+1]>DZERO) ||
1114 (DEQUAL(spar1[1]+tref2,sval2[1]+tref2) && sp1inf[ki+1]<DZERO) ||
1115 (DEQUAL(spar1[0]+tref1,sval1[0]+tref1) && sp1inf[ki]>DZERO) ||
1116 (DEQUAL(spar1[0]+tref1,sval1[1]+tref1) && sp1inf[ki]<DZERO)
1117
1118 )
1119 kstat = 1;
1120 }
1121 krem1 = kstat;
1122
1123
1124
1125 /* Check if curve between start and midd point cross the
1126 boundary */
1127
1128 s1305(spar1,siparmid,sval1,sval2,&kbound,sipar1,&kstat);
1129 if (kstat<0) goto error;
1130 krem2 = kstat;
1131
1132 /* We now have the following cases:
1133 kstat == 0 : Line between epar1 and epar2 outside,
1134 If this happens when kdir=1, then
1135 just forget the start point. If it happens
1136 when kdir=2, then we just stop the marching.
1137 kstat == 1 : Line between epar1 and epar2 inside.
1138 Continue iteration.
1139 kstat == 2 : We step out of the patch. Clip to the edge
1140 of the patch. Update start point.
1141 kstat == 3 : We step into the patch. Clip to the edge
1142 of the patch. Update endpoint
1143 kstat == 4 : We go from the boundary and out. Try next
1144 iteration direction.
1145
1146 */
1147 if (krem1==0 || krem2==0)
1148 {
1149 if (kdir==1) knbinf--;
1150 goto nextdir;
1151 }
1152 else if ((krem1 !=1 || krem2 !=1) && krem1 != 4 && krem2 !=4)
1153 {
1154
1155 /* If we clip to the boundary, forget any guide point
1156 identified. The actual action is dependent on which
1157 of the points spar2 or sipar1 indicates the crossing */
1158
1159 kstat1 = 0;
1160 if (krem2==2 || krem2==3)
1161 {
1162 s9clipimp(spar1,siparmid,ps1,eimpli,ideg,sval1,sval2,
1163 aepsge,sipnt1,sipar1,&kstat);
1164 if (kstat<0) goto error;
1165 if (krem2==3 && kstat==1) kstpch = 4;
1166 kstat1 = kstat;
1167 krem = krem2;
1168 }
1169 if (kstat1 != 1 && (krem1 ==2 || krem1==3))
1170 {
1171 s9clipimp(spar1,spar2,ps1,eimpli,ideg,sval1,sval2,
1172 aepsge,sipnt1,sipar1,&kstat);
1173 if (kstat<0) goto error;
1174 if (krem1==3 && kstat==1) kstpch = 4;
1175 kstat1 = kstat;
1176 krem = krem1;
1177 }
1178
1179 if (kstat1 == 1)
1180 {
1181 /* Check that the relationship between the two surfaces
1182 * has not been interchanged, by making the cross product
1183 * of the normal vectors in current point and the point
1184 * found by iteration. Then make the scalar product of
1185 * these vectors. If the scalar product is negative then
1186 * we have either jumped to another branch or passed a
1187 * singularity, iterprete this as the iteration has diverged
1188 * In addition we don't want the direction of the tangents
1189 * change to much. We set a limit of approximately 41 degrees
1190 * by testing on a cosin value of 0.75
1191 * Make normal vectors in implicit surface for both points
1192 * Make also sure that the curve in the parameter plane
1193 * does not turn more than 90 degrees.
1194 */
1195
1196 s1331(spnt1,eimpli,ideg,-1,tval,snorm1,&kstat);
1197 if (kstat < 0) goto error;
1198 s1331(sipnt1,eimpli,ideg,-1,tval,snorm2,&kstat);
1199 if (kstat < 0) goto error;
1200
1201 if(ideg == 1003 || ideg == 1004 || ideg == 1005)
1202 {
1203 tdum = s6scpr(snorm1,snorm2,kdim);
1204 }
1205 else
1206 {
1207 s6crss(spnt1+ksizem3,snorm1,sdum1);
1208 (void)s6norm(sdum1,kdim,sdum1,&kstat);
1209 if (kstat < 0) goto error;
1210 s6crss(sipnt1+ksizem3,snorm2,sdum2);
1211 (void)s6norm(sdum2,kdim,sdum2,&kstat);
1212 if (kstat < 0) goto error;
1213
1214 tdum = s6scpr(sdum1,sdum2,kdim);
1215 }
1216
1217
1218 /* Check that sipar1 lies on the same side of spar1 as
1219 the tangent at spar1 */
1220
1221 s6diff(sipar1,spar1,2,sdum1);
1222 tdump = s6scpr(sdum1,sp1inf+7*(knbinf-1)+2,2);
1223 /* The test below was added to detect the case when the
1224 normal of the parametric surface turns */
1225 if(s6scpr(spnt1+ksizem3,sipnt1+ksizem3,kdim) < 0.0)
1226 {
1227 tdum = -tdum;
1228 }
1229 }
1230
1231 /* An intersection point has only been found when kstat==1
1232 */
1233 if (kstat1==1 && tdump >= (double)0.0 && tdum > (double)0.75)
1234 {
1235 /* If krem=3 we step into the patch, if krem=2 we step
1236 out of the patch */
1237
1238 if (krem==2 || krem==3)
1239 {
1240 /* Since we clip, set kstpch=0, no guide point reached */
1241
1242 /* If krem==3 we step into the patch, make new
1243 start point of segment */
1244
1245 if (krem==3) knbinf--;
1246
1247 memcopy(spnt2,sipnt1,ksize,DOUBLE);
1248 memcopy(spar2,sipar1,2,DOUBLE);
1249
1250 s1306(sipnt1,sipar1,eimpli,ideg,s3dinf+10*knbinf,
1251 sp1inf+7*knbinf,&kstat);
1252 if (kstat<0) goto error;
1253
1254
1255 /* Make sure that the tangents of previous and the new point
1256 point in the same direction */
1257
1258 if (knbinf>0)
1259 tdum = s6scpr(s3dinf+10*(knbinf-1)+3,
1260 s3dinf+10*knbinf+3,kdim);
1261 else
1262 tdum = s6scpr(startg,s3dinf+3,kdim);
1263
1264
1265 if (tdum < DZERO)
1266 {
1267 /* Change tangent direction 3-D and in
1268 parameter plane
1269 */
1270 sp1 = s3dinf + 10*knbinf + 3;
1271 sp1[0] = -sp1[0];
1272 sp1[1] = -sp1[1];
1273 sp1[2] = -sp1[2];
1274 sp1 = sp1inf + 7*knbinf + 2;
1275 sp1[0] = -sp1[0];
1276 sp1[1] = -sp1[1];
1277 }
1278
1279 /* If the new end point tangent points out go
1280 to next direction */
1281
1282 ki = 7*knbinf;
1283 if( (sp1inf[ki+1] <= sval2[0] && sp1inf[ki+3] < DZERO) ||
1284 (sp1inf[ki+1] >= sval2[1] && sp1inf[ki+3] > DZERO) ||
1285 (sp1inf[ki ] <= sval1[0] && sp1inf[ki+2] < DZERO) ||
1286 (sp1inf[ki ] >= sval1[1] && sp1inf[ki+2] > DZERO) )
1287 {
1288 knbinf++;
1289 goto nextdir;
1290 }
1291 else if(krem == 2 &&
1292 ((sp1inf[ki+1] <= sval2[0] && sp1inf[ki+3] >= DZERO) ||
1293 (sp1inf[ki+1] >= sval2[1] && sp1inf[ki+3] <= DZERO) ||
1294 (sp1inf[ki ] <= sval1[0] && sp1inf[ki+2] >= DZERO) ||
1295 (sp1inf[ki ] >= sval1[1] && sp1inf[ki+2] <=DZERO) ))
1296 {
1297 /* We were marching out ou the patch, but the tangent
1298 points in, half step length */
1299 kstpch = 0;
1300 }
1301
1302
1303
1304 }
1305 }
1306
1307 else
1308 {
1309 /* Divergence or point on wrong side in the parameter
1310 plane or 3-d */
1311 kstpch = 0;
1312 koutside_resolution = 0;
1313 }
1314 }
1315 else if (krem1==4 || krem2==4)
1316 goto nextdir;
1317
1318 }
1319
1320 /* Update step length if new endpoint is to be found */
1321
1322 if (kstpch==0)
1323 {
1324 if (tdist<DZERO)
1325 {
1326 tnew = tstep/(double)10.0;
1327 }
1328 else
1329 {
1330 tfak = MAX(tdist/aepsge,(double)1.0);
1331 tfak = (double)2.0*pow(tfak,ONE_FOURTH);
1332 tnew = MIN(tstep/(double)2.0,tstep/tfak);
1333 }
1334 if (DEQUAL(tmax+tnew,tmax+tstep)) goto nextdir;
1335 tstep = tnew;
1336 }
1337 }
1338
1339 /* If kstpch= -1,1,3 or 4 then a point is accepted and
1340 * snxt1 points to the position and derivatives
1341 * of the accepted point. */
1342
1343 /* If we have accepted a segment pointing in the opposite
1344 * direction of the previous segment, something very wrong
1345 * has happened, and we go out with an error. */
1346 if (knbinf >= 2)
1347 {
1348 s6diff(s3dinf + 10*knbinf, s3dinf + 10*(knbinf - 1), kdim, sdiffcur);
1349 s6diff(s3dinf + 10*(knbinf-1), s3dinf + 10*(knbinf - 2), kdim, sdiffprev);
1350 if (s6scpr(sdiffcur, sdiffprev ,kdim) < DZERO)
1351 {
1352 /* We have a problem with degeneracy, quit now. */
1353 goto war03;
1354 }
1355 /* printf("%7.13f\n", s6scpr(sdiffcur, sdiffprev ,kdim)); */
1356 }
1357
1358 /* Update number of intersection points */
1359
1360 knbinf++;
1361
1362 /* Copy point and parameter pair descriptions */
1363
1364 memcopy(spnt1,spnt2,ksize,DOUBLE);
1365 memcopy(spar1,spar2,2,DOUBLE);
1366
1367 /* Update guide point pointers */
1368
1369 if (kstpch == 1)
1370 {
1371 kguide++;
1372 kgdir = 1;
1373
1374 /* Test if end of guide point array reached */
1375
1376 if (kguide >= kpoint) goto nextdir;
1377
1378 }
1379 if (kstpch == -1)
1380 {
1381 kguide--;
1382 kgdir = -1;
1383
1384 /* Test if start of guide point array reached */
1385
1386 if (1 >= kguide) goto nextdir;
1387 }
1388
1389 /* Make new radius of curvature */
1390
1391 trad = *(s3dinf + 10*knbinf - 1);
1392 tstep = s1311(trad,aepsge,tmax,&kstat);
1393 if (kstat<0) goto error;
1394
1395 /* Test if start point reached, e.g. that the curve is closed */
1396
1397 if (kstpch == 3)
1398 {
1399 /* Closed curve found */
1400 goto finished;
1401 }
1402
1403
1404 /* End while loop */
1405 }
1406 /* End inside */
1407 /* INSIDE TEST REMOVED BECAUSE OF CHANGED STRATEGY
1408 }
1409 */
1410 nextdir:;
1411 /* End two step directions */
1412 }
1413
1414 finished:
1415
1416 /* In certain cases too many marched point may be found. These cases are:
1417
1418 - Open curve and start of marching first guide point
1419 - Open curve and start of marching last guide point
1420 - Closed curve and this found in second marching direction
1421
1422 In these cases some of the found points have to be discarded */
1423
1424 scorpnt = s3dinf;
1425 scorpar = sp1inf;
1426
1427 if (kstpch !=3 && kpoint>1)
1428 {
1429
1430 /* Open curve */
1431
1432 if ( (kstart==1 && kgd1 == kpoint) ||
1433 (kstart==kpoint && kgd1==1) )
1434 {
1435 /* First marching direction traced curve */
1436
1437 knbinf = knb1;
1438 }
1439 else if ( (kstart==1 && kguide==kpoint) ||
1440 (kstart==kpoint && kguide==1) )
1441 {
1442 /* Second marching direction traced curve */
1443
1444 scorpnt = scorpnt + 10*(knb1-1);
1445 scorpar = scorpar + 7*(knb1-1);
1446 knbinf = knbinf - knb1 + 1;
1447 }
1448 }
1449 else if (kpoint>1)
1450 {
1451 /* Closed curve, correct if result of second marching direction */
1452
1453 if (kdir != 1)
1454 {
1455 /* Second marching direction, disc ard result of first direction */
1456
1457 scorpnt = scorpnt + 10*(knb1-1);
1458 scorpar = scorpar + 7*(knb1-1);
1459 knbinf = knbinf - knb1 + 1;
1460 }
1461 }
1462
1463 /* A curve is traced out only if at least two points are found, if less
1464 points found try to pick out constant parameter line */
1465
1466 interpolate:
1467
1468 if (knbinf>1)
1469 {
1470 if (igraph == 1 && knbinf > 1)
1471 {
1472 /* Output curve through s6line and s6move */
1473
1474 s6move(scorpnt);
1475 for (ki=1,sp1=scorpnt+10;ki<knbinf;ki++,sp1+=10)
1476 s6line(sp1);
1477 }
1478
1479 if (icur > 0 && knbinf >1)
1480 {
1481
1482 /* Make 3-D representation of intersection curve */
1483
1484 kpar = 0;
1485
1486 /* We allocate space for parametrization array */
1487
1488
1489 spar = newarray(knbinf,DOUBLE);
1490 if (spar == SISL_NULL) goto err101;
1491
1492 s1359(scorpnt,aepsge,kdim,knbinf,kpar,spar,&q3dcur,&kstat);
1493 if (kstat < 0) goto error;
1494
1495 /* Set pointer in intcurve object to 3-D curve */
1496 pintcr -> pgeom = q3dcur;
1497
1498 if (icur == 2)
1499 {
1500 /* Make curve in parameter plane */
1501
1502 kdim = 2;
1503 kpar = 1;
1504 s1359(scorpar,aepsge,kdim,knbinf,kpar,spar,&qp1cur,&kstat);
1505 if (kstat < 0) goto error;
1506
1507
1508 /* Set pointersin intcurve object to curves in parameter plane */
1509 if (kpar1 ==2)
1510 {
1511 pintcr -> ppar1 = qp1cur;
1512 }
1513 else
1514 {
1515 pintcr -> ppar2 = qp1cur;
1516 }
1517 }
1518 }
1519 }
1520 /* Dont use s9constline for the silhouette curves -- it won't work! */
1521 else if(pintcr->ipoint > 1 && ideg < 1003)
1522 {
1523
1524 /* If no points produced on intersection curve */
1525
1526 s1313_s9constline(ps1,eimpli,ideg,aepsge,pintcr,
1527 icur,igraph,&kstat);
1528 if (kstat<0) goto error;
1529 if (kstat==0) goto err185;
1530 }
1531 else
1532 goto err185;
1533
1534 if (kdiv==1) goto war03;
1535 *jstat = 0;
1536
1537 goto out;
1538
1539 /* Iteration can not continue */
1540 war03: *jstat = 3;
1541 goto out;
1542
1543 /* Error in space allocation */
1544 err101: *jstat = -101;
1545 s6err("s1313",*jstat,kpos);
1546 goto out;
1547
1548 /* Error in surface description parameter direction does not exist */
1549 err123: *jstat = -123;
1550 s6err("s1313",*jstat,kpos);
1551 goto out;
1552
1553
1554 /* Error - SISL_NULL pointer was given */
1555 err150 :
1556 *jstat = -150;
1557 s6err("s1313",*jstat,kpos);
1558 goto out;
1559
1560
1561 /* Only degenerate or singular guide points */
1562 err185: *jstat = -185;
1563 goto out;
1564
1565
1566 /* Error in lower leve function */
1567 error:
1568 *jstat = kstat;
1569 s6err("s1313",*jstat,kpos);
1570 goto out;
1571
1572 out:
1573
1574 /* Free allocated space */
1575
1576 if (sgpar != SISL_NULL) freearray(sgpar);
1577 if (sgd1 != SISL_NULL) freearray(sgd1);
1578 if (s3dinf != SISL_NULL) freearray(s3dinf);
1579 if (sp1inf != SISL_NULL) freearray(sp1inf);
1580 if (spar != SISL_NULL) freearray(spar);
1581
1582
1583 return;
1584 }
1585
1586 #if defined(SISLNEEDPROTOTYPES)
1587 static void
s1313_s9constline(SISLSurf * ps1,double eimpli[],int ideg,double aepsge,SISLIntcurve * pintcr,int icur,int igraph,int * jstat)1588 s1313_s9constline(SISLSurf *ps1,double eimpli[],int ideg,
1589 double aepsge,SISLIntcurve *pintcr,int icur,
1590 int igraph,int *jstat)
1591 #else
1592 static void s1313_s9constline(ps1,eimpli,ideg,aepsge,pintcr,
1593 icur,igraph,jstat)
1594 SISLSurf *ps1;
1595 double eimpli[];
1596 int ideg;
1597 double aepsge;
1598 SISLIntcurve *pintcr;
1599 int icur;
1600 int igraph;
1601 int *jstat;
1602 #endif
1603 /*
1604 *********************************************************************
1605 *
1606 * PURPOSE : To check if the parameter pairs describe an intersection
1607 * curve that is a constant parameter line in the parameter
1608 * plane of a surface and to produce the description of
1609 * the curve according to the specifications.
1610 * DO NOT USE IT FOR SILHOUETTE CURVES -- IT WON'T WORK!
1611 *
1612 *
1613 * INPUT : ps1 - Pointer to surface.
1614 * eimpli - Description of the implicit surface
1615 * ideg - The degree of the implicit surface
1616 * ideg=1: Plane
1617 * ideg=2; Quadric surface
1618 * ideg=1001: Torus surface
1619 * ideg=1003: Silhouette line parallel projection
1620 * ideg=1004: Silhouette line perspective projection
1621 * ideg=1005: Silhouette line circular projection
1622 * aepsco - Computational resolution.
1623 * aepsge - Geometry resolution.
1624 * amax - Maximal allowed step length. Not used.
1625 * icur - Indicator telling if a 3-D curve is to be made
1626 * 0 - Don't make 3-D curve
1627 * 1 - Make 3-D curve
1628 * 2 - Make 3-D curve and curves in parameter plane
1629 * igraph - Indicator telling if the curve is to be outputted
1630 * through function calls:
1631 * 0 - don't output curve through function call
1632 * 1 - output as straight line segments through
1633 * s6move and s6line
1634 *
1635 *
1636 *
1637 * INPUT/OUTPUT:pintcr - The intersection curve. When comming as input
1638 * only parameter values in the parameter plane
1639 * exist. When comming as output the 3-D geometry
1640 * and possibly the curve in the parameter plane
1641 * of the surface are added.
1642 *
1643 * OUTPUT: jstat - status messages
1644 * = 1 : Constant parameter line is intersection
1645 * = 0 : No intersection along constant parameter
1646 * line.
1647 * < 0 : error
1648 *
1649 *
1650 * METHOD :
1651 * REFERENCES :
1652 *
1653 *
1654 *-
1655 * CALLS :
1656 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway, 2. July 1989
1657 * Revised by : Mike Floater, SI, 1991-01
1658 * Tried to add perspective and circular silhouettes (ideg=1004,ideg=1005)
1659 * but more work is needed. After improving s1313, s1309, and s1331
1660 * for all silhouette
1661 * types -- ideg=1003,1004,1005 this routine is no longer
1662 * useable even for ideg=1003. The problem is that there is not
1663 * enough information in qc2 and its first derivative for the
1664 * calls to s1331 and s1309. Perhaps the solution is to write a
1665 * new routine specially for silhouette curve intersections.
1666 * Until this problem is fixed S1313_S9CONSTLINE MUST NOT BE
1667 * CALLED WHEN ideg=1003, ideg=1004, ideg=1005 (see s1313).
1668 *
1669 *********************************************************************
1670 */
1671 {
1672 int ki,kj,kl; /* Control variables in for loops */
1673 int kk,kn,kk1,kn1,kk2,kn2;/* Orders and numbers of knots */
1674 int kpoint; /* Number of points in guide curve */
1675 int kleft = 0; /* Pointer into knot vector */
1676 int kpar1; /* Number of parameter direction in 1st. obj */
1677 int kpar2; /* Number of parameter direction in 2st. obj */
1678 int ktype; /* Type of intersection curve */
1679 int kpos=0; /* Position of error */
1680 int kstat; /* Status variable returned form routine */
1681 int kdir=0; /* constant parameter line direction */
1682 int knbpnt; /* Number of points on constant parameter line */
1683 double tmax1,tmin1; /* Minimum and maximum of first comp of guide points */
1684 double tmax2,tmin2; /* Minimum and maximum of first comp of guide points */
1685 double tmax; /* Maximum 3-D SISLbox side */
1686 double tdist,tang; /* Distance and angle error */
1687 double *st,*st1,*st2; /* Pointers to knot vectors */
1688 double *spoint; /* Pointer to points on constant parameter line */
1689 double *sp1; /* Pointer into array */
1690 double tsize1,tsize2; /* Length of knot intervals */
1691 double sval1[2]; /* Limits of parameter plane in first SISLdir */
1692 double sval2[2]; /* Limits of parameter plane in second SISLdir */
1693 double sder[6]; /* SISLPoint and derivative on curve */
1694 double sider[3]; /* SISLPoint on implicit surface */
1695 double snorm[3]; /* Normal on implicit surface */
1696 double tsumold,tsum,tval;/* Parameter values */
1697 double *sgpar1=SISL_NULL; /* Parameter pairs of guide point in surf 1 */
1698 double *sgpar2=SISL_NULL; /* Parameter pairs of guide point in surf 2 */
1699 SISLCurve *qc1=SISL_NULL; /* Pointer to 3-D curve */
1700 SISLCurve *qc2=SISL_NULL; /* Pointer to 3-D curve */
1701
1702 SISLCurve *qp1cur=SISL_NULL; /* Pointer to curve in first parameter plane*/
1703
1704
1705
1706 /* Make maximal step length based on box-size of surface */
1707
1708 sh1992su(ps1,0,aepsge,&kstat);
1709 if (kstat < 0) goto error;
1710
1711 tmax = MAX(ps1->pbox->e2max[0][0] - ps1->pbox->e2min[0][0],
1712 ps1->pbox->e2max[0][1] - ps1->pbox->e2min[0][1]);
1713 tmax = MAX(tmax,ps1->pbox->e2max[0][2] - ps1->pbox->e2min[0][2]);
1714
1715 /* Find a none singular start point for the marching process */
1716
1717 kpoint = pintcr->ipoint;
1718 kpar1 = pintcr->ipar1;
1719 kpar2 = pintcr->ipar2;
1720 sgpar1 = pintcr->epar1;
1721 sgpar2 = pintcr->epar2;
1722 ktype = pintcr->itype;
1723
1724
1725 /* Initiate pointers to intersection curve and intersection curve in
1726 parameter plane */
1727
1728 pintcr -> pgeom = SISL_NULL;
1729 pintcr -> ppar1 = SISL_NULL;
1730 pintcr -> ppar2 = SISL_NULL;
1731
1732
1733 /* Initiate parameter direction boundaries */
1734 kk1 = ps1 -> ik1;
1735 kn1 = ps1 -> in1;
1736 st1 = ps1 -> et1;
1737 sval1[0] = st1[kk1-1];
1738 sval1[1] = st1[kn1];
1739 kk2 = ps1 -> ik2;
1740 kn2 = ps1 -> in2;
1741 st2 = ps1 -> et2;
1742 sval2[0] = st2[kk2-1];
1743 sval2[1] = st2[kn2];
1744
1745
1746 /* Test that first object has 2 parameter
1747 direction and second object 0 */
1748
1749 if (kpar1 == 2 && kpar2 == 0)
1750 {
1751 /* Everithing is ok */
1752 ;
1753 }
1754 else if (kpar1 == 0 && kpar2 == 2)
1755 {
1756 sgpar1 = sgpar2;
1757 }
1758 else
1759 {
1760 goto err123;
1761 }
1762
1763
1764 /* Run through the parameter pairs to decide if a constant parameter line
1765 is possible */
1766
1767 tmax1 = tmin1 = sgpar1[0];
1768 tmax2 = tmin2 = sgpar1[1];
1769
1770 for (ki=1,kj=2,kl=3 ; ki < kpoint ; ki++,kj+=2,kl+=2)
1771 {
1772 tmin1 = MIN(tmin1,sgpar1[kj]);
1773 tmax1 = MAX(tmax1,sgpar1[kj]);
1774 tmin2 = MIN(tmin2,sgpar1[kl]);
1775 tmax2 = MAX(tmax2,sgpar1[kl]);
1776 }
1777
1778 tsize1 = st1[kn1] - st1[kk1-1];
1779 tsize2 = st2[kn2] - st2[kk2-1];
1780
1781 /* Check if constant parameter value within tolerance */
1782
1783 if (DEQUAL((tmin1+tsize1),(tmax1+tsize1)) )
1784 {
1785 /* Intersection possible constant parameter line with first parameter
1786 constant value constant.
1787
1788 1. Pick out curve from surface
1789 2. Pick out relevant part of curve */
1790
1791 kdir = 1;
1792
1793 s1437(ps1,(tmin1+tmax1)/2.0,&qc1,&kstat);
1794 if (kstat < 0) goto error;
1795
1796 s1712(qc1,tmin2,tmax2,&qc2,&kstat);
1797 if (kstat < 0) goto error;
1798 }
1799 else if (DEQUAL((tmin2+tsize2),(tmax2+tsize2)) )
1800 {
1801 /* Intersection possible constant parameter line with first parameter
1802 constant value constant.
1803
1804 1. Pick out curve from surface
1805 2. Pick out relevant part of curve */
1806
1807 kdir = 2;
1808
1809 s1436(ps1,(tmin2+tmax2)/2.0,&qc1,&kstat);
1810 if (kstat < 0) goto error;
1811
1812 s1712(qc1,tmin1,tmax1,&qc2,&kstat);
1813 if (kstat < 0) goto error;
1814 }
1815 else
1816 goto war00;
1817
1818 st = qc2 -> et;
1819 kk = qc2 -> ik;
1820 kn = qc2 -> in;
1821
1822 /* Run through 2*kn points of the curve and check that they lie in the
1823 implicit surface by calculating the 2*kn points. */
1824
1825 tsumold = st[kk-1];
1826
1827 for (ki=0 ; ki <kn ; ki++)
1828 {
1829 if (kk>1)
1830 {
1831 /* Make parameter value to use for calculation of curve point */
1832
1833 for (kl=1,kj=ki+1,tsum=(double)0.0 ; kl<kk ; kl++)
1834 tsum += st[kj++];
1835
1836 tsum = tsum/(kk-1);
1837 }
1838 else
1839 tsum = st[ki];
1840
1841 tval = tsum;
1842
1843 for (kj=0 ; kj<2 ; kj++)
1844 {
1845 /* Calculate point on curve */
1846
1847 s1221(qc2,1,tval,&kleft,sder,&kstat);
1848 if (kstat < 0) goto error;
1849
1850 /* Calculate normal to implicit surface */
1851
1852 s1331(sder,eimpli,ideg,-1,sider,snorm,&kstat);
1853 if (kstat < 0) goto error;
1854
1855 /* Project point onto implicit surface */
1856
1857 tdist = fabs(s1309(sder,snorm,eimpli,ideg,&kstat));
1858 if (kstat<0) goto error;
1859 if (kstat==2) goto war00;
1860
1861 /* Both the position of the two points should be within the relative
1862 computer resolution for the point to be accepted. Correspondingly the
1863 direction of the intersection curve and the constant parameter line
1864 should be within the computer resolution to be accepted. */
1865
1866 if (DNEQUAL(tdist+tmax,tmax))
1867 goto war00;
1868
1869 /* Distance within tolerance, check that the angle between surface
1870 normal and curve tangent is PIHALF, if both these vectors have a
1871 nonzero length. */
1872
1873 if (s6length(snorm,3,&kstat) != (double)0.0 &&
1874 s6length(sder+3,3,&kstat) != (double)0.0 )
1875 {
1876 tang = s6ang(snorm,sder+3,3);
1877 if (DNEQUAL(fabs(tang),PIHALF) ) goto war00;
1878 }
1879 tval = (tsumold+tsum)/(double)2.0;
1880 }
1881 tsumold = tsum;
1882 }
1883
1884 /* Intersection curve along constant parameter line, make right actions
1885 concerning drawing and/or creation of the curve */
1886
1887 if (igraph == 1)
1888 {
1889 /* Draw curve, first break into straight line segments */
1890
1891 s1605(qc2,aepsge,&spoint,&knbpnt,&kstat);
1892 if (kstat < 0) goto error;
1893
1894 if (knbpnt>1)
1895 {
1896 /* Draw curve */
1897
1898 s6move(spoint);
1899 for (ki=1,sp1=spoint+3 ; ki<knbpnt ; ki++,sp1+=3)
1900 s6line(sp1);
1901 }
1902 freearray(spoint);
1903 }
1904
1905 if (icur >= 1)
1906 {
1907 /* Set pointer to 3-D curve */
1908
1909 pintcr -> pgeom = qc2;
1910 qc2 = SISL_NULL;
1911 }
1912
1913 if (icur == 2)
1914 {
1915 /* Make curve in parameter plane */
1916
1917 double svert[4],sknot[4];
1918
1919 if (kdir==1)
1920 {
1921 svert[0] = svert[2] = (tmin1+tmax1)/(double)2.0;
1922 svert[1] = tmin2;
1923 svert[3] = tmax2;
1924 sknot[0] = sknot[1] = tmin2;
1925 sknot[2] = sknot[3] = tmax2;
1926 }
1927 else
1928 {
1929 svert[0] = tmin1;
1930 svert[2] = tmax1;
1931 svert[1] = svert[3] = (tmin2+tmax2)/(double)2.0;
1932 sknot[0] = sknot[1] = tmin1;
1933 sknot[2] = sknot[3] = tmax1;
1934 }
1935 qp1cur = newCurve(2,2,sknot,svert,1,2,1);
1936 if (qp1cur==SISL_NULL) goto err101;
1937 pintcr -> ppar1 = qp1cur;
1938 }
1939
1940 /* war01: */
1941 *jstat = 1;
1942 goto out;
1943
1944 /* Iteration can not continue */
1945 war00: *jstat = 0;
1946 goto out;
1947
1948
1949 /* Error in space allocation */
1950 err101: *jstat = -101;
1951 s6err("s1313",*jstat,kpos);
1952 goto out;
1953
1954 /* Error in surface description parameter direction does not exist */
1955 err123: *jstat = -123;
1956 s6err("s1313_s9constline",*jstat,kpos);
1957 goto out;
1958
1959 /* Error in lower leve function */
1960 error:
1961 *jstat = kstat;
1962 s6err("s1313_s9constline",*jstat,kpos);
1963 goto out;
1964
1965 out:;
1966 if (qc1 != SISL_NULL) freeCurve(qc1);
1967 if (qc2 != SISL_NULL) freeCurve(qc2);
1968 }
1969