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 #define S1912
43 
44 #include "sislP.h"
45 
46 #if defined(SISLNEEDPROTOTYPES)
47 typedef void (*fparamProc)(double[], int[], int, int, int, double, double *,
48 			   double *[], double *[], int *);
49 typedef void (*fknotsProc)(double[], int, int, int, double *[], int *);
50 #else
51 typedef void (*fparamProc)();
52 typedef void (*fknotsProc)();
53 #endif
54 
55 
56 #if defined(SISLNEEDPROTOTYPES)
57 void
s1912(fparamProc fparam,fknotsProc fknots,double econd[],int ntype[],int inpt,double astpar,int ik,int idim,int iopen,double * cendpar,SISLCurve ** rcurve,double ** gpar,int * jnbpar,int * jstat)58    s1912 (fparamProc fparam, fknotsProc fknots,
59        double econd[], int ntype[], int inpt, double astpar, int ik,
60        int idim, int iopen, double *cendpar, SISLCurve ** rcurve,
61        double **gpar, int *jnbpar, int *jstat)
62 #else
63 void
64    s1912 (fparam, fknots, econd, ntype, inpt, astpar, ik, idim, iopen,
65        cendpar, rcurve, gpar, jnbpar, jstat)
66      fparamProc fparam;
67      fknotsProc fknots;
68      double econd[];
69      int ntype[];
70      int inpt;
71      double astpar;
72      int ik;
73      int idim;
74      int iopen;
75      double *cendpar;
76      SISLCurve **rcurve;
77      double **gpar;
78      int *jnbpar;
79      int *jstat;
80 #endif
81 /*
82 *********************************************************************
83 *
84 * PURPOSE    : Compute a B-spline curve interpolating a set of points.
85 *              The points can be assigned derivative conditions. The
86 *              curve can be open, closed, or closed and periodic.
87 *
88 * INPUT      : fparam - Routine computing parametrization of point set.
89 *	       fknots - Routine coputing knot vector of point set.
90 *              econd  - Array of interpolation conditions. Dimension
91 *                       is inpt*idim.
92 *              ntype  - Array containing kind of condition. Dimension
93 *                       is inpt.
94 *                       =  0 : A point is given.
95 *                       =  d : The d'th derivatative condition to the
96 *                              previous point is given.
97 *                       = -d : The d'th derivatative condition to the
98 *                              next point is given.
99 *              inpt   - Number of interpolation conditions.
100 *              astpar - Start parameter of parametrization.
101 *              ik     - Order of interpolating curve.
102 *              idim   - Dimension of geometry space.
103 *              iopen - Indicates if the curve is to be open, closed or
104 *                       periodic.
105 *
106 * OUTPUT     : cendpar - End parameter of parametrization.
107 *              rcurve  - Interpolating curve.
108 *	       gpar    - The distinct parameter values.
109 *	       jnbpar  - Number of distinct parameter values.
110 *              jstat   - status messages
111 *                        = 1      : Specified parametrization method
112 *                                   replaced by cord length parametrization.
113 *                                         = 0      : ok
114 *
115 * METHOD     :
116 *
117 * REFERENCES :
118 *
119 * CALLS      :	s1905, s1891, s1713, s1750
120 *
121 * WRITTEN BY : Vibeke Skytt, SI, 91-04.
122 * REVISED BY : Trond Vidar Stensby, SI, 91-07
123 * REWISED BY : Vibeke Skytt, 94-03. Changed the concept of closed,
124 *                                   non-periodic.
125 * REWISED BY : Johannes Kaasa, 95-11. Fixed error in output of the parameter
126 *              values (included the last parameter value).
127 *
128 *********************************************************************
129 */
130 {
131   int kstat = 0;		/* Status variable.                             */
132   int kpos = 0;
133   int ki;			/* Counter.                                     */
134   int knpt;			/* Number of accepted interpolation conditions. */
135   int kn;			/* Number of coefficients of B-spline curve.    */
136   int kordr;			/* Local order of curve.                        */
137   int kright = 1;		/* One equation system to solve in interpolation. */
138   int knlr = 0;			/* Indicates shape of interpolation matrix.     */
139   int knrc = 0;			/* Indicates shape of interpolation matrix.     */
140   int kopen;                    /* Local open/closed parameter. Closed,
141 				   non-periodic is treated as an open curve.*/
142   int *ltype = SISL_NULL;		/* Type of accepted interpolation conditions.   */
143   double *scond = SISL_NULL;		/* Array containing interpolation conditions.   */
144   double *spar1 = SISL_NULL;		/* Parametrization array of interpolation conditions. */
145   double *spar2 = SISL_NULL;		/* Parametrization array used to make knot vector. */
146   double *sknot = SISL_NULL;		/* Knot vector of curve.                           */
147   double *scoef = SISL_NULL;		/* Coefficients of curve.                          */
148   int *sder = SISL_NULL;		/* Vector of derivative indicators.                */
149   SISLCurve *qc = SISL_NULL;		/* Interpolation curve.                            */
150   SISLCurve *qc2 = SISL_NULL;	/* Interpolation curve.                            */
151 
152   *jstat = 0;
153 
154   /* Set local open/closed parameter. */
155 
156   kopen = (iopen == SISL_CRV_PERIODIC) ? 0 : 1;
157 
158   /* Test interpolation conditions, and adjust the input conditions
159      if necessary.  */
160 
161   s1905 (econd, ntype, inpt, ik, idim, iopen, &scond, &ltype, &knpt, &kstat);
162   if (kstat < 0)
163     goto error;
164 
165   /* Set local order.  */
166 
167   kordr = MIN (ik, knpt);
168 
169   /* Allocate scratch for derivative indicator. */
170 
171   if ((sder = newarray (knpt, INT)) == SISL_NULL)
172     goto err101;
173 
174   for (ki = 0; ki < knpt; ki++)
175     sder[ki] = (int) fabs ((double) ltype[ki]);
176 
177   /* Compute parametrization of point set. */
178 
179   (* fparam) (scond, ltype, knpt, idim, kopen, astpar, cendpar,
180 	      &spar1, &spar2, &kstat);
181   if (kstat < 0) goto error;
182 
183   /* Make knot vector of curve.  */
184 
185   if (iopen == SISL_CRV_PERIODIC)
186     {
187       /* Closed, periodic curve. */
188 
189       knlr = kordr / 2;
190       knrc = kordr - knlr - 1;
191       knpt--;
192     }
193 
194   /* Produce knot vector. */
195 
196   (* fknots) (spar1, knpt, kordr, kopen, &sknot, &kstat);
197   if (kstat < 0) goto error;
198 
199   /* Perform interpolation.  */
200 
201   s1891 (spar1, scond, idim, knpt, kright, sder, kopen, sknot,
202 	 &scoef, &kn, kordr, knlr, knrc, &kstat);
203   if (kstat < 0) goto error;
204 
205   /* Express the curve as a curve object.  */
206 
207   qc = newCurve (kn, kordr, sknot, scoef, 1, idim, 1);
208   if (qc == SISL_NULL) goto err101;
209 
210   qc->cuopen = iopen;
211 
212   if (kordr < ik)
213     {
214       /* The order of the curve is less than expected. Increase the order. */
215 
216       qc2 = SISL_NULL;
217       s1750 (qc, ik, &qc2, &kstat);
218       if (kstat < 0) goto error;
219 
220       if (qc != SISL_NULL) freeCurve (qc);
221       qc = qc2;
222     }
223 
224   /* Interpolation performed. */
225 
226   /* Find distinct parameter values. */
227 
228   *gpar = spar1;
229   *jnbpar = 0;
230 
231   for (ki = 1; spar1[ki] < *cendpar; ki++)
232     {
233       if (spar1[ki - 1] < spar1[ki])
234 	(*gpar)[(*jnbpar)++] = spar1[ki-1];
235     }
236   (*gpar)[(*jnbpar)++] = spar1[ki-1];
237   (*gpar)[(*jnbpar)++] = spar1[ki];
238 
239   *gpar = increasearray (*gpar, *jnbpar, DOUBLE);
240 
241   *rcurve = qc;
242   goto out;
243 
244 
245   /* Error in scratch allocation.  */
246 
247 err101:
248   *jstat = -101;
249   s6err ("s1912", *jstat, kpos);
250   goto out;
251 
252   /* Error in lower level routine. */
253 
254 error:
255   *jstat = kstat;
256   s6err ("s1912", *jstat, kpos);
257   goto out;
258 
259 out:
260   /* Free scratch occupied by local arrays. */
261 
262   if (spar2 != SISL_NULL)
263     freearray (spar2);
264   if (scond != SISL_NULL)
265     freearray (scond);
266   if (scoef != SISL_NULL)
267     freearray (scoef);
268   if (sknot != SISL_NULL)
269     freearray (sknot);
270   if (sder != SISL_NULL)
271     freearray (sder);
272   if (ltype != SISL_NULL)
273     freearray (ltype);
274 
275   return;
276 }
277