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: s1955.c,v 1.3 2001-03-19 15:58:57 afr Exp $
45 *
46 */
47
48
49 #define S1955
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1955(SISLCurve * pc1,SISLCurve * pc2,double aepsco,double aepsge,int * jpt,double ** gpar1,double ** gpar2,int * jcrv,SISLIntcurve *** wcurve,int * jstat)55 s1955(SISLCurve *pc1,SISLCurve *pc2,double aepsco,double aepsge,
56 int *jpt,double **gpar1,double **gpar2,int *jcrv,SISLIntcurve ***wcurve,int *jstat)
57 #else
58 void s1955(pc1,pc2,aepsco,aepsge,jpt,gpar1,gpar2,jcrv,wcurve,jstat)
59 SISLCurve *pc1;
60 SISLCurve *pc2;
61 double aepsco;
62 double aepsge;
63 int *jpt;
64 double **gpar1;
65 double **gpar2;
66 int *jcrv;
67 SISLIntcurve ***wcurve;
68 int *jstat;
69 #endif
70 /*
71 *********************************************************************
72 *
73 *********************************************************************
74 *
75 * PURPOSE : Find the closest point between the two curve pc1 and pc2.
76 *
77 *
78 *
79 * INPUT : pc1 - Pointer to the first curve in the closest point
80 * problem.
81 * pc2 - Pointer to the second curve in the closest point
82 * problem.
83 * aepsco - Computational resolution.
84 * aepsge - Geometry resolution.
85 *
86 *
87 *
88 * OUTPUT : jpt - Number of single closest points.
89 * gpar1 - Array containing the parameter values of the
90 * single closest points in the parameter intervals
91 * of the first curve. The points lie continuous.
92 * Closest curves are stored in wcurve.
93 * gpar2 - Array containing the parameter values of the
94 * single closest points in the parameter intervals
95 * of the second curve. The points lie continuous.
96 * jcrv - Number of closest curves.
97 * wcurve - Array containing descriptions of the closest
98 * curves. The curves are only described by points
99 * in the parameter area. The curve-pointers points
100 * to nothing. (See descrjption of Intcurve
101 * in intcurve.dcl).
102 * If the curves given as input are degnenerate an
103 * intersection point can be returned as an intersection
104 * curve. Use s1327 to decide if an intersection curve
105 * is a point on one of the curves.
106 * jstat - status messages
107 * > 0 : warning
108 * = 0 : ok
109 * < 0 : error
110 *
111 *
112 * METHOD : Generate the difference surface between the two curves,
113 * and find those points of this surface closest to origo.
114 *
115 *
116 *
117 * REFERENCES :
118 *
119 *-
120 * CALLS : s1954 - Find closest points between a surface and a point.
121 * s1956 - Create difference surface between two curves.
122 * freeSurf - Free space occupied by a surface.
123 * newIntcurve - Create new closest curve.
124 *
125 * WRITTEN BY : Vibeke Skytt, SI, 88-11.
126 * Revised by : Paal Fugelli, SINTEF, Oslo, Norway, 1994-08. Updated to handle
127 * input curves with periodic basis correctly.
128 *
129 *********************************************************************
130 */
131 {
132 int kstat = 0; /* Local status variable. */
133 int kpos = 0; /* Position of error. */
134 int ki,kj; /* Counters. */
135 int kdim; /* Dimension of the space in which the
136 curves lie. */
137 double *sorigo = SISL_NULL; /* Array representing origo. */
138 double *spar1 = SISL_NULL; /* Parameter-values in the parameter
139 interval of the first curve of the
140 closest curve. */
141 double *spar2 = SISL_NULL; /* Parameter-values in the parameter
142 interval of the second curve of the
143 closest curve. */
144 double *ssing = SISL_NULL; /* Single closest points between
145 difference surface and origo. */
146 SISLSurf *qsdiff = SISL_NULL; /* Difference surface. */
147 SISLIntcurve *qcurve; /* Pointer to closest curve. */
148 SISLCurve *qkreg1 = SISL_NULL; /* First incrv with ensured k-regular basis. */
149 SISLCurve *qkreg2 = SISL_NULL; /* Second incrv with ensured k-regular basis. */
150 int k1=1,k2=2,k4=4; /* Constants */
151
152
153 /* Ensure k-regular basis for input curves. */
154
155 if ( pc1 -> cuopen == SISL_CRV_PERIODIC )
156 {
157 /* Cyclic (i.e. periodic) basis. */
158
159 make_cv_kreg(pc1, &qkreg1, &kstat);
160 if ( kstat < 0 ) goto error;
161 }
162 else
163 qkreg1 = pc1;
164
165 if ( pc2 -> cuopen == SISL_CRV_PERIODIC )
166 {
167 /* Cyclic (i.e. periodic) basis. */
168
169 make_cv_kreg(pc2, &qkreg2, &kstat);
170 if ( kstat < 0 ) goto error;
171 }
172 else
173 qkreg2 = pc2;
174
175
176 /* Test input. */
177
178 kdim = qkreg1 -> idim;
179 if ( kdim != qkreg2 -> idim ) goto err106;
180
181 /* Allocate space for the point origo. */
182
183 if ( (sorigo = new0array(kdim, DOUBLE)) == SISL_NULL ) goto err101;
184
185 /* Generate the difference surface between the curves. */
186
187 s1956(qkreg1, qkreg2, &qsdiff, &kstat);
188 if ( kstat < 0 ) goto error;
189
190 if ( kstat > 0 )
191 {
192
193 /* The curves define a closest interval. Allocate space
194 for points defining the interval. */
195
196 if ( (spar1 = newarray(k2, DOUBLE)) == SISL_NULL ) goto err101;
197 if ( (spar2 = newarray(k2, DOUBLE)) == SISL_NULL ) goto err101;
198
199 if ( kstat == 1 )
200 {
201
202 /* The curves are equal. */
203
204 spar1[0] = *(qkreg1->et + qkreg1->ik - 1);
205 spar1[1] = *(qkreg1->et + qkreg1->in);
206 spar2[0] = *(qkreg2->et + qkreg2->ik - 1);
207 spar2[1] = *(qkreg2->et + qkreg2->in);
208
209 }
210 else if ( kstat == 2 )
211 {
212
213 /* The curves are equal, but oppositely directed. */
214
215 spar1[0] = *(qkreg1->et + qkreg1->ik - 1);
216 spar1[1] = *(qkreg1->et + qkreg1->in);
217 spar2[0] = *(qkreg2->et + qkreg2->in);
218 spar2[1] = *(qkreg2->et + qkreg2->ik - 1);
219 }
220
221 *wcurve = SISL_NULL;
222 if ( (*wcurve = newarray(k1, SISLIntcurve*)) == SISL_NULL ) goto err101;
223 *jcrv = 1;
224
225 **wcurve = SISL_NULL;
226 if ( (**wcurve = newIntcurve(k2,k1,k1,spar1,spar2,k4)) == SISL_NULL ) goto err101;
227
228 }
229 else
230 {
231
232 /* Find the closest points between the difference function and origo.*/
233
234 s1954(qsdiff, sorigo, kdim, aepsco, aepsge, jpt, &ssing, jcrv, wcurve, &kstat);
235 if ( kstat < 0 ) goto error;
236
237 /* Handle periodicity (remove extraneous points) */
238
239 if ( (*jpt) > 1 && kdim > 1 && (pc1 -> cuopen == SISL_CRV_PERIODIC ||
240 pc2 -> cuopen == SISL_CRV_PERIODIC) )
241 {
242 for ( ki=0; ki < (*jpt); ki++ )
243 {
244 if ( (pc1 -> cuopen == SISL_CRV_PERIODIC &&
245 ssing[2*ki] == pc1 -> et[pc1->in]) ||
246 (pc2 -> cuopen == SISL_CRV_PERIODIC &&
247 ssing[2*ki+1] == pc2 -> et[pc2->in]) )
248 {
249 (*jpt)--;
250 ssing[2*ki] = ssing[2*(*jpt)];
251 ssing[2*ki+1] = ssing[2*(*jpt)+1];
252 ki--;
253 }
254 }
255 }
256
257 if ( (*jpt) > 0 )
258 {
259
260 /* Allocate space for output of single closest points. */
261
262 *gpar1 = *gpar2 = SISL_NULL;
263 if ( (*gpar1 = newarray(*jpt, DOUBLE)) == SISL_NULL ) goto err101;
264 if ( (*gpar2 = newarray(*jpt, DOUBLE)) == SISL_NULL ) goto err101;
265
266 /* Copy single closest points to output arrays. */
267
268 for ( ki=0; ki < (*jpt); ki++ )
269 {
270 *(*gpar1 + ki) = ssing[2*ki];
271 *(*gpar2 + ki) = ssing[2*ki+1];
272 }
273 }
274
275 for ( kj=0; kj < (*jcrv); kj++ )
276 {
277
278 /* Correct parameter-arrays of closest curve. */
279
280 qcurve = *(*wcurve + kj);
281 qcurve -> ipar1 = qcurve -> ipar2 = 1;
282
283 /* Allocate space for new parameter arrays. */
284
285 spar1 = spar2 = SISL_NULL;
286 if ( (spar1 = newarray(qcurve->ipoint, DOUBLE)) == SISL_NULL ) goto err101;
287 if ( (spar2 = newarray(qcurve->ipoint, DOUBLE)) == SISL_NULL ) goto err101;
288
289 /* Copy parameter-values to new arrays. */
290
291 for ( ki=0; ki < qcurve->ipoint; ki++ )
292 {
293 spar1[ki] = *(qcurve->epar1+2*ki);
294 spar2[ki] = *(qcurve->epar1+2*ki+1);
295 }
296
297 freearray(qcurve->epar1);
298 qcurve -> epar1 = spar1;
299 qcurve -> epar2 = spar2;
300 }
301 }
302
303 /* Closest points found. */
304
305 *jstat = 0;
306 goto out;
307
308
309 /* Error in space allocation. */
310
311 err101:
312 *jstat = -101;
313 s6err("s1955",*jstat,kpos);
314 goto out;
315
316 /* Error in input. Dimensions conflicing. */
317
318 err106:
319 *jstat = -106;
320 s6err("s1955",*jstat,kpos);
321 goto out;
322
323 /* Error in lower level routine. */
324
325 error:
326 *jstat = kstat;
327 s6err("s1955",*jstat,kpos);
328 goto out;
329
330 out:
331
332 /* Free space occupied local variables. */
333
334 if ( qkreg1 && qkreg1 != pc1 ) freeCurve(qkreg1);
335 if ( qkreg2 && qkreg2 != pc2 ) freeCurve(qkreg2);
336 if (ssing) freearray(ssing);
337 if (sorigo) freearray(sorigo);
338 if (qsdiff) freeSurf(qsdiff);
339
340 return;
341 }
342