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: s1379.c,v 1.3 2001-03-19 15:58:48 afr Exp $ 45 * 46 */ 47 48 49 #define S1379 50 51 #include "sislP.h" 52 53 #if defined(SISLNEEDPROTOTYPES) 54 void s1379(double ep[],double ev[],double epar[],int im,int idim,SISLCurve ** rcurve,int * jstat)55 s1379(double ep[],double ev[],double epar[],int im,int idim, 56 SISLCurve **rcurve,int *jstat) 57 #else 58 void s1379(ep,ev,epar,im,idim,rcurve,jstat) 59 double ep[]; 60 double ev[]; 61 double epar[]; 62 int im; 63 int idim; 64 SISLCurve **rcurve; 65 int *jstat; 66 #endif 67 /* 68 ************************************************************************ 69 * 70 * Purpose: To compute the cubic Hermit interpolant to the data given 71 * by the points ep, the derivatives ev and paramterization epar. 72 * The curve is represented as a B-spline curve.If the first and 73 * last points are exactly equal (down to the last bit) the a periodic 74 * basis with the first 1 single knot, then a 3 tupple knot is made 75 * at the start and 3 tupple knot followed by a single knot at the 76 * end. If also the derivatives at the start and end are equal then 77 * all knots will be double. 78 * 79 * 80 * Input: 81 * ep - Array containing the point in sequence 82 * (x,y,..,x,y,..) 83 * ev - Array containing the derivatives in sequence 84 * (x,y,..,x,y,..) 85 * epar - Parametrization array. The array should be increasing 86 * in value. 87 * im - Number of point and derivatives 88 * idim - The dimension of the space the points and derivatives 89 * lie in 90 * Output: 91 * rcurve - Pointer to the curve produced 92 * jstat - Status variable 93 * < 0 - Error. 94 * Method: 95 * The knot vector will have 4-tupple, 3-tupple or 2-tupple knots at 96 * epar[0] and epar[im-1]. This is decided by checking if the input 97 * data is cyclic: 98 * if first point != last point 4-tupple knots 99 * if first == last and there derivatives different 3-tupple 100 * if both position and derivatives equal 2-tupple 101 * In the case not 4-tupple knots the new knot intervals are made cyclic. 102 * 103 * Suppose we have reached data point no. j which corresponds to the 104 * parameter value z=epar(j), i.e., the knot vector et has a double 105 * knot at z, and these knots must be et(2*j+1) and et(2*j+2). 106 * then there are only two B-splines which are nonzero at z, namely 107 * B(2*j-1) and B(2*j) (remember everything is cubic), and these are 108 * also the only B-splines with nonzero derivative at z. we can 109 * therefore determine the coefficients of these two B-splines, 110 * ec(2*j-1) and ec(2*j), by solving a 2x2 linear system OF 111 * equations, and this can be done directly. 112 * Suppose the distance from z to the previous knot is h1 (measured 113 * in parameter interval, so h1= et(2*j+2)-et(2*j), and the distance 114 * to the next knot is h2 = et(2*j+3) -et(2*j+1). then the two 115 * B-splines and their derivatives have the following values at z: 116 * 117 * B(2*j-1,z)=h2/(h1+h2), B(2*j,z)=h1/(h1+h2) 118 * dB(2*j-1,z)=-3/(h1+h2), dB(2*j,z)=3/(h1+h2). 119 * 120 * solving the linear system for ec(2*j-1) and ec(2*j) yields 121 * 122 * ec(2*j-1)=ep(j) - h1*ev(j)/3, 123 * ec(2*j)=ep(j) + h2*ev(j)/3. 124 * 125 * it can also be easily checked that these equations are valid for 126 * the first two and last two coefficients as well, provided one 127 * sets h1=0 and h2=0, respectively. 128 * 129 * REFERENCES : 130 * 131 *- 132 * CALLS : 133 * 134 * WRITTEN BY : Tor Dokken, SI 1988-11 135 * Revised by : Tor Dokken, SI 1992-04 136 * Revised by : Johs. Kaasa,SI 1993-01 137 * 138 ********************************************************************* 139 */ 140 { 141 int ki,kj; /* Loop variables */ 142 int kk; /* Polynomial order */ 143 int kn; /* Number of vertices */ 144 int kpoint; /* Pointer into point and derivative array */ 145 int kcoef; /* Pointer into coefficient array */ 146 int kpos=0; /* Position of error */ 147 int kthis; /* Current point */ 148 int kstat=0; /* Status variable */ 149 int kcycpos = 1; /* Flag telling if first and last points are equal */ 150 int kcycder = 1; /* Flag telling if first and last derviatives are equal */ 151 double *st=SISL_NULL; /* Knot vector */ 152 double *scoef=SISL_NULL; /* B-spline vertices */ 153 double th1,th2; /* Parameter intervals */ 154 155 156 157 /* Check input */ 158 159 if (im < 2) goto err181; 160 if (idim < 1) goto err102; 161 162 /* Set the dimension and order of the spline space */ 163 164 kn = 2*im; 165 kk = 4; 166 167 /* Allocate arrays for temporary storage of knots and vertices */ 168 169 st = newarray(kn+kk,DOUBLE); 170 if (st == SISL_NULL) goto err101; 171 scoef = newarray(idim*kn,DOUBLE); 172 if (scoef == SISL_NULL) goto err101; 173 174 /* Check if the curve is periodic, e.g. if first and last points are 175 equal and/or that first and last derivates are equal */ 176 177 /* for (kj=0, kcycpos=1 ; kj<idim && kcycpos == 1 ; kj++) 178 if (ep[kj] != ep[idim*(im-1)+kj]) kcycpos =0; */ 179 for (kj=0, kcycpos=1 ; kj<idim && kcycpos == 1 ; kj++) 180 if (DNEQUAL(ep[kj], ep[idim*(im-1)+kj])) kcycpos =0; 181 182 /* for (kj=0, kcycder=1 ; kj<idim && kcycder == 1 ; kj++) 183 if (ev[kj] != ev[idim*(im-1)+kj]) kcycder= 0; */ 184 for (kj=0, kcycder=1 ; kj<idim && kcycder == 1 ; kj++) 185 if (DNEQUAL(ev[kj], ev[idim*(im-1)+kj])) kcycder= 0; 186 187 /* Make the knot vector, first all knots except the two first and the two last */ 188 189 for (ki=2,kj=0 ; ki<kn+2 ; ki+=2, kj++) 190 st[ki] = st[ki+1] = epar[kj]; 191 192 193 194 /* Make the two first and two last knots */ 195 196 if (kcycder == 1 && kcycpos == 1) 197 { 198 /* Two first knots to be shifted */ 199 200 st[0]= st[1] = epar[0] - (epar[im-1]-epar[im-2]); 201 st[kn+2]= st[kn+3] = epar[im-1] + epar[1] - epar[0]; 202 } 203 else if (kcycder ==0 && kcycpos ==1) 204 { 205 /* First and last knot to be shifted */ 206 207 st[0] = epar[0] - (epar[im-1]-epar[im-2]); 208 st[1] = st[2]; 209 st[kn+2] = st[kn]; 210 st[kn+3] = epar[im-1] + epar[1] - epar[0]; 211 } 212 else 213 { 214 /* k-regular basis */ 215 216 st[0] = st[1] = st[2]; 217 st[kn+2] = st[kn+3] = st[kn]; 218 } 219 220 /* Compute knot vector and coefficients as indicated above */ 221 222 for (kj=0, kcoef=0, kpoint = 0 ; kj < kn ; kj+=2) 223 { 224 th1 = st[kj+3] - st[kj+1]; 225 th2 = st[kj+4] - st[kj+2]; 226 227 /* Compute coefficient no kj */ 228 229 kthis = kpoint; 230 for (ki=0;ki<idim;ki++,kpoint++) 231 { 232 scoef[kcoef++] = ep[kpoint] - th1*ONE_THIRD*ev[kpoint]; 233 } 234 235 /* Compute coefficient no kj+1 */ 236 237 kpoint = kthis; 238 for (ki=0;ki<idim;ki++,kpoint++) 239 { 240 scoef[kcoef++] = ep[kpoint] + th2*ONE_THIRD*ev[kpoint]; 241 } 242 } 243 244 /* Make new curve object */ 245 246 *rcurve = newCurve(kn,kk,st,scoef,1,idim,1); 247 if (*rcurve == SISL_NULL) goto err101; 248 249 /* Remove unneccesarry knots */ 250 251 s6crvcheck(*rcurve,&kstat); 252 if (kstat<0) goto error; 253 254 /* Periodicity flag */ 255 if (kcycpos) 256 { 257 test_cyclic_knots((*rcurve)->et,(*rcurve)->in,(*rcurve)->ik,&kstat); 258 if (kstat<0) goto error; 259 if (kstat == 2) (*rcurve)->cuopen = SISL_CRV_PERIODIC; 260 } 261 262 /* Calculation completed */ 263 264 *jstat = 0; 265 goto out; 266 267 268 /* Error in space allocation. Return zero. */ 269 270 271 /* Error in space allocation */ 272 err101: *jstat = -101; 273 s6err("s1379",*jstat,kpos); 274 goto out; 275 276 277 /* Dimension less than 1*/ 278 err102: *jstat = -102; 279 s6err("s1379",*jstat,kpos); 280 goto out; 281 282 /* Too few interpolation conditions */ 283 284 err181: *jstat = -181; 285 s6err("s1379",*jstat,kpos); 286 goto out; 287 288 error: *jstat =kstat; 289 s6err("s1379",*jstat,kpos); 290 goto out; 291 292 out: 293 if (st != SISL_NULL) freearray(st); 294 if (scoef != SISL_NULL) freearray(scoef); 295 296 return; 297 } 298