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 #define S1357 44 45 #include "sislP.h" 46 47 #if defined(SISLNEEDPROTOTYPES) 48 void s1357(double epoint[],int inbpnt,int idim,int ntype[],double epar[],int icnsta,int icnend,int iopen,int ik,double astpar,double * cendpar,SISLCurve ** rc,double ** gpar,int * jnbpar,int * jstat)49 s1357(double epoint[],int inbpnt,int idim,int ntype[],double epar[], 50 int icnsta,int icnend,int iopen,int ik,double astpar, 51 double *cendpar,SISLCurve **rc,double **gpar,int *jnbpar,int *jstat) 52 #else 53 void s1357(epoint,inbpnt,idim,ntype,epar,icnsta,icnend,iopen,ik,astpar, 54 cendpar,rc,gpar,jnbpar,jstat) 55 double epoint[]; 56 int inbpnt; 57 int idim; 58 int ntype[]; 59 double epar[]; 60 int icnsta; 61 int icnend; 62 int iopen; 63 int ik; 64 double astpar; 65 double *cendpar; 66 SISLCurve **rc; 67 double **gpar; 68 int *jnbpar; 69 int *jstat; 70 #endif 71 /* 72 ************************************************************************* 73 * 74 * Purpose: To calculate a B-spline curve interpolating a set of points. 75 * The points can be assigned a tangent (derivative). 76 * The curve can be closed or open. If end-conditions are 77 * conflicting, the condition closed curve rules out other 78 * end conditions. 79 * The parametrization is given by the array Epar. 80 * 81 * Input: 82 * Epoint - Array (length idim*inbpnt) containing the points/ 83 * derivatives to be interpolated. 84 * Inbpnt - No. of points/derivatives in the epoint array. 85 * Idim - The dimension of the space in which the points lie. 86 * ntype - Array (length inbpnt) containing type indicator for 87 * points/derivatives/second-derivatives: 88 * 1 - Ordinary point. 89 * 2 - Knuckle point. (Is treated as an ordinary point.) 90 * 3 - Derivative to next point. 91 * 4 - Derivative to prior point. 92 * ( 5 - Second derivative to next point. ) 93 * ( 6 - Second derivative to prior point. ) 94 * 13 - Start-point of tangent to next point. 95 * 14 - End-point of tangent to prior point. 96 * Epar - Array containing the wanted parametrization. Only parameter 97 * values corresponding to position points are given. 98 * For closed curves, one additional parameter value 99 * must be spesified. The last entry contains 100 * the parametrization of the repeted start point. 101 * (if the endpoint is equal to the startpoint of 102 * the interpolation the lenght of the array should 103 * be equal to inpt1 also in the closed case). 104 * Icnsta - Additional condition at the start of the curve: 105 * 0 : No additional condition. 106 * 1 : Zero curvature at start. 107 * Icnend - Additional condition at the end of the curve: 108 * 0 : No additional condition. 109 * 1 : Zero curvature at end. 110 * Iopen - Flag telling if the curve should be open or closed: 111 * 1 : The curve should be open. 112 * 0 : The curve should be closed. 113 * -1 : The curve should be closed and periodic. 114 * Ik - The order of the B-spline curve to be produced. 115 * Astpar - Parameter-value to be used at the start of the curve. 116 * 117 * Output: 118 * Jstat - status variable: 119 * < 0 : Error. 120 * = 0 : Ok. 121 * > 0 : Warning. 122 * cendpar - Parameter-value used at the end of the curve. 123 * Rc - Pointer to output-curve. 124 * Gpar - Pointer to the parameter-values of the points in 125 * the curve. Represented only once, although derivatives 126 * and second-derivatives will have the same parameter- 127 * value as the points. 128 * Jnbpar - No. of different parameter-values. 129 * 130 * Method: First the parametrization of the curve and the type 131 * specification of points/derivatives are checked and/or 132 * corrected. Then the knots are calculated, and the 133 * interpolation is performed. 134 * 135 * Calls: s1907,s1908,s1902,s1891,s1713,s1750,s6err. 136 * 137 * Written by: Trond Vidar Stensby, SI, 1991-07 138 * The Fortran routine, p19539, is written by: T. Dokken SI. 139 * Bug fix: Michael Floater, 82.9.93. The construction of gpar and *jnbpar 140 * at the end of the routine was incorrect. 141 * REWISED BY: Vibeke Skytt, 03.94. This routine corresponds to s1358, 142 * but differ in the use of the parameter 143 * iopen and the input array ntype is of 144 * type int. 145 * REWISED BY: Johannes Kaasa, 01.98. Made sure the start and end points 146 * are correct in the periodic case. I also made proper 147 * handling when the order is higher than the number of 148 * interpolation points. 149 ***************************************************************** 150 */ 151 { 152 int kpos = 0; 153 SISLCurve *qc = SISL_NULL; /* Temporary SISLCurves. */ 154 SISLCurve *qc2 = SISL_NULL; 155 SISLCurve *dummy = SISL_NULL; 156 int *ltype = SISL_NULL; /* Type of interpolation condition 157 (temporary) */ 158 int *ltype2 = SISL_NULL; /* Type of interpolation condition (finial) */ 159 int *sder = SISL_NULL; /* Vector of derivative indicators. */ 160 int knpt; /* Number of interpolation conditions. */ 161 int kordr; /* Local order. */ 162 int kstat; /* Status variable. */ 163 int kn; /* Number of coefficients of B-spline curve.*/ 164 int ki, kj, kl; /* Loop control variable. */ 165 int kright = 1; /* One equation system to solve in 166 interpolation. */ 167 int knlr = 0; /* Indicates shape of interpolation matrix. */ 168 int knrc = 0; /* Indicates shape of interpolation matrix. */ 169 int kopen; /* Local open/closed parameter. Closed, 170 non-periodic is treated as an open curve.*/ 171 int kpair; /* Pair order or not. */ 172 int kcont; /* Continuity in start/end point. */ 173 int kleft; /* Pointer into knot array. */ 174 int klast1, klast2; /* Last element in array. */ 175 double split_par; /* Splitting parameter. */ 176 double tdiff; /* Length of parameter interval. */ 177 double *lpar = SISL_NULL; /* Parameter values. (temporary) */ 178 double *lcond = SISL_NULL; /* Interpolation conditions. (temporary) */ 179 double *sknot = SISL_NULL; /* Knot vector. */ 180 double *spar = SISL_NULL; /* Parameter valued. (finial) */ 181 double *scond = SISL_NULL; /* Interpolation conditions. (finial) */ 182 double *scoef = SISL_NULL; /* Coefficients of curve. */ 183 double *temp = SISL_NULL; /* Temporary storage. */ 184 185 *jstat = 0; 186 187 /* If necessary reduce the order. */ 188 189 kordr = MIN (ik, inbpnt); 190 191 /* Set local open/closed parameter. */ 192 193 kopen = (iopen == SISL_CRV_PERIODIC) ? 0 : 1; 194 195 /* Check pair order or not. */ 196 197 kpair = (kordr % 2 == 0) ? 1 : 0; 198 199 /* Transform interpolation conditions. */ 200 201 s1907 (epoint, ntype, epar, iopen, icnsta, icnend, inbpnt, 202 idim, &lcond, <ype, &lpar,&knpt, &kstat); 203 if (kstat < 0) goto error; 204 205 /* Test interpolation conditions, and adjust the input conditions 206 if necessary. */ 207 208 s1908 (lcond, ltype, lpar, knpt, kordr, idim, iopen, &scond, <ype2, 209 &spar, &knpt, &kstat); 210 if (kstat < 0) goto error; 211 212 /* Allocate scratch for derivative indicator. */ 213 214 sder = newarray (knpt, INT); 215 if (sder == SISL_NULL) goto err101; 216 217 for (ki = 0; ki < knpt; ki++) 218 sder[ki] = abs (ltype2[ki]); 219 220 221 if (iopen == SISL_CRV_PERIODIC) 222 { 223 knlr = kordr / 2; 224 knrc = kordr - knlr - 1; 225 knpt--; 226 } 227 228 /* Produce knot vector. */ 229 230 s1902 (lpar, knpt, kordr, kopen, &sknot, &kstat); 231 if (kstat < 0) goto error; 232 233 /* Perform interpolation. */ 234 235 s1891 (spar, scond, idim, knpt, kright, sder, kopen, sknot, 236 &scoef, &kn, kordr, knlr, knrc, &kstat); 237 if (kstat < 0) goto error; 238 239 /* Corrected start and end for odd order periodic curves. */ 240 241 if (iopen == SISL_CRV_PERIODIC && kpair) 242 { 243 if ((temp = newarray(idim, double)) == SISL_NULL) 244 goto err101; 245 246 /* Find number of present starting knots. */ 247 248 ki = 0; 249 while (sknot[ki] < epar[0]) ki++; 250 251 klast1 = kn + kordr - 1; 252 klast2 = idim*(kn - 1); 253 tdiff = lpar[knpt] - lpar[0]; 254 for ( ; ki < (kordr - 1); ki++) 255 { 256 temp[0] = sknot[knpt - 1] - tdiff; 257 for (kj = klast1; kj > 0; kj--) 258 sknot[kj] = sknot[kj - 1]; 259 sknot[0] = temp[0]; 260 261 for (kl = 0; kl < idim; kl++) 262 temp[kl] = scoef[(knpt - 1)*idim + kl]; 263 for (kj = klast2; kj > 0; kj -= idim) 264 { 265 for (kl = 0; kl < idim; kl++) 266 scoef[kj + kl] = scoef[kj - idim + kl]; 267 } 268 for (kl = 0; kl < idim; kl++) 269 scoef[kl] = temp[kl]; 270 } 271 } 272 273 /* Express the curve as a curve object. */ 274 275 qc = newCurve (kn, kordr, sknot, scoef, 1, idim, 1); 276 if (qc == SISL_NULL) goto err101; 277 qc->cuopen = iopen; 278 279 /* Corrected start and end for even order periodic curves. */ 280 281 if (iopen == SISL_CRV_PERIODIC && !kpair && 282 fabs(qc->et[qc->ik - 1] - epar[0]) > REL_PAR_RES) 283 { 284 285 /* Shift the start/end point. */ 286 287 split_par = epar[0]; 288 while (split_par < qc->et[qc->ik-1]) 289 split_par += (qc->et[qc->in] - qc->et[qc->ik-1]); 290 while (split_par > qc->et[qc->in]) 291 split_par -= (qc->et[qc->in] - qc->et[qc->ik-1]); 292 s1710(qc, split_par, &qc2, &dummy, &kstat); 293 if (kstat != 2) goto error; 294 tdiff = qc2->et[qc2->ik - 1] - epar[0]; 295 if (fabs(tdiff) > REL_PAR_RES) 296 { 297 for (ki = 0; ki < (qc2->in + qc2->ik); ki++) 298 qc2->et[ki] -= tdiff; 299 } 300 301 /* qc2 is represented on a closed basis, and flagged 302 as SISL_CRV_CLOSED. We have to open it again. First 303 find the continuity in the new end point. */ 304 305 kcont = qc->ik - 1 - max(s6knotmult(qc->et, qc->ik, qc->in, 306 &kleft, epar[0], &kstat), 1); 307 if (kstat < 0) goto error; 308 309 make_cv_cyclic(qc2, kcont, &kstat); 310 if (kstat < 0) goto error; 311 312 if (qc != SISL_NULL) freeCurve (qc); 313 qc = qc2; 314 qc2 = SISL_NULL; 315 316 } 317 318 if (kordr < ik) 319 { 320 /* The order of the curve is less than expected. Increase the order. */ 321 322 qc2 = SISL_NULL; 323 s1750 (qc, ik, &qc2, &kstat); 324 if (kstat < 0) goto error; 325 326 if (qc != SISL_NULL) freeCurve (qc); 327 qc = qc2; 328 qc2 = SISL_NULL; 329 } 330 331 /* Set open/closed parameter of curve. */ 332 333 qc->cuopen = iopen; 334 335 /* Set end of parameter interval. */ 336 337 *cendpar = *(qc->et + qc->in); 338 339 /* Interpolation performed. */ 340 341 /* Find distinct parameter values. */ 342 343 *gpar = lpar; 344 345 /* Bug fix, MSF, 28.9.93. Change 346 *jnbpar = 1; 347 to: */ 348 349 *jnbpar = 0; 350 351 for (ki = 1; ki<knpt; ki++) 352 { 353 if (spar[ki - 1] < spar[ki]) 354 355 /* Bug fix. MSF, 28.9.93 Change (*gpar)[(*jnbpar)++] = spar[ki]; to: */ 356 357 (*gpar)[(*jnbpar)++] = spar[ki-1]; 358 } 359 360 /* Bug fix. MSF, 28.9.93 Change (*gpar)[(*jnbpar)++] = spar[ki]; to: */ 361 362 (*gpar)[(*jnbpar)++] = spar[ki-1]; 363 364 *rc = qc; 365 goto out; 366 367 368 /* Error in scratch allocation. */ 369 370 err101: 371 *jstat = -101; 372 s6err ("s1357", *jstat, kpos); 373 goto out; 374 375 /* Error in lower level routine. */ 376 377 error: 378 *jstat = kstat; 379 s6err ("s1357", *jstat, kpos); 380 goto out; 381 382 out: 383 /* Free scratch occupied by local arrays. */ 384 385 if (scond != SISL_NULL) freearray (scond); 386 if (scoef != SISL_NULL) freearray (scoef); 387 if (sder != SISL_NULL) freearray (sder); 388 if (ltype != SISL_NULL) freearray (ltype); 389 if (ltype2 != SISL_NULL) freearray (ltype2); 390 if (lcond != SISL_NULL) freearray (lcond); 391 if (sknot != SISL_NULL) freearray (sknot); 392 if (spar != SISL_NULL) freearray (spar); 393 if (temp != SISL_NULL) freearray (temp); 394 395 return; 396 } 397