//=========================================================================== // GoTools Core - SINTEF Geometry Tools Core library, version 2.0.1 // // Copyright (C) 2000-2007, 2010 SINTEF ICT, Applied Mathematics, Norway. // // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License // as published by the Free Software Foundation version 2 of the License. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software // Foundation, Inc., // 59 Temple Place - Suite 330, // Boston, MA 02111-1307, USA. // // Contact information: E-mail: tor.dokken@sintef.no // SINTEF ICT, Department of Applied Mathematics, // P.O. Box 124 Blindern, // 0314 Oslo, Norway. // // Other licenses are also available for this software, notably licenses // for: // - Building commercial software. // - Building software whose source code you wish to keep private. //=========================================================================== 00018 #ifndef __BORLANDC__ 00019 #include "boost/lambda/lambda.hpp" 00020 #endif 00021 00022 using namespace std; 00023 #ifndef __BORLANDC__ 00024 using namespace boost::lambda; 00025 #endif 00026 00027 namespace Go 00028 { 00029 #ifdef __BORLANDC__ 00030 namespace 00031 { 00033 class ScaleBy : public std::unary_function<double, double> 00034 { 00035 double m_scale; 00036 00037 public: 00038 ScaleBy(const double& scale) : m_scale(scale) {} 00039 00040 double operator()(const double& value) { return m_scale * value; } 00041 }; 00042 } // anonymous namespace 00043 #endif 00044 00045 //=========================================================================== 00046 // this is an attempt to optimize the point evaluation routine. 00047 // If it is not satisfactory for some reason, the original routine exists commented-out 00048 // below 00049 void SplineSurface::point(Point& result, double upar, double vpar) const 00050 //=========================================================================== 00051 { 00052 result.resize(dim_); 00053 const int uorder = order_u(); 00054 const int vorder = order_v(); 00055 const int unum = numCoefs_u(); 00056 int kdim = rational_ ? dim_ + 1 : dim_; 00057 00058 static ScratchVect<double, 10> Bu(uorder); 00059 static ScratchVect<double, 10> Bv(vorder); 00060 static ScratchVect<double, 4> tempPt(kdim); 00061 static ScratchVect<double, 4> tempResult(kdim); 00062 00063 Bu.resize(uorder); 00064 Bv.resize(vorder); 00065 tempPt.resize(kdim); 00066 tempResult.resize(kdim); 00067 00068 // compute tbe basis values and get some data about the spline spaces 00069 basis_u_.computeBasisValues(upar, Bu.begin()); 00070 basis_v_.computeBasisValues(vpar, Bv.begin()); 00071 const int uleft = basis_u_.lastKnotInterval(); 00072 const int vleft = basis_v_.lastKnotInterval(); 00073 00074 // compute the tensor product value 00075 const int start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim; 00076 00077 register double* ptemp; 00078 register const double* co_ptr = rational_ ? &rcoefs_[start_ix] : &coefs_[start_ix]; 00079 fill(tempResult.begin(), tempResult.end(), double(0)); 00080 00081 for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) { 00082 register const double bval_v = *bval_v_ptr; 00083 fill(tempPt.begin(), tempPt.end(), 0); 00084 for (register double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) { 00085 register const double bval_u = *bval_u_ptr; 00086 for (ptemp = tempPt.begin(); ptemp != tempPt.end(); ++ptemp) { 00087 *ptemp += bval_u * (*co_ptr++); 00088 } 00089 } 00090 ptemp = tempPt.begin(); 00091 for (register double* p = tempResult.begin(); p != tempResult.end(); ++p) { 00092 *p += (*ptemp++) * bval_v; 00093 } 00094 co_ptr += kdim * (unum - uorder); 00095 } 00096 00097 copy(tempResult.begin(), tempResult.begin() + dim_, result.begin()); 00098 if (rational_) { 00099 const double w_inv = double(1) / tempResult[kdim - 1]; 00100 #ifdef __BORLANDC__ // C++Builder does not support boost lambda 00101 transform(result.begin(), result.end(), result.begin(), ScaleBy(w_inv)); 00102 #else 00103 transform(result.begin(), result.end(), result.begin(), _1 * w_inv); 00104 #endif 00105 } 00106 } 00107 00108 // this was the original point evaluation routine. 00109 00110 // //=========================================================================== 00111 // void SplineSurface::point(Point& result, double upar, double vpar) const 00112 // //=========================================================================== 00113 // { 00114 // if (result.dimension() != dim_) 00115 // result.resize(dim_); 00116 00117 // // Take care of the rational case 00118 // const std::vector<double>& co = rational_ ? rcoefs_ : coefs_; 00119 // int kdim = dim_ + (rational_ ? 1 : 0); 00120 00121 // // Make temporary storage for the basis values and a temporary 00122 // // computation cache. 00123 // Go::ScratchVect<double, 10> b0(basis_u_.order()); 00124 // Go::ScratchVect<double, 10> b1(basis_v_.order()); 00125 // Go::ScratchVect<double, 4> temp(kdim); 00126 // Go::ScratchVect<double, 4> restemp(kdim); 00127 // std::fill(restemp.begin(), restemp.end(), 0.0); 00128 00129 // // Compute the basis values and get some data about the spline spaces 00130 // basis_u_.computeBasisValues(upar, &b0[0]); 00131 // int uleft = basis_u_.lastKnotInterval(); 00132 // int uorder = basis_u_.order(); 00133 // int unum = basis_u_.numCoefs(); 00134 // basis_v_.computeBasisValues(vpar, &b1[0]); 00135 // int vleft = basis_v_.lastKnotInterval(); 00136 // int vorder = basis_v_.order(); 00137 00138 // // Compute the tensor product value 00139 // int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); 00140 // for (int jj = 0; jj < vorder; ++jj) { 00141 // std::fill(temp.begin(), temp.end(), 0.0); 00142 // for (int ii = 0; ii < uorder; ++ii) { 00143 // for (int dd = 0; dd < kdim; ++dd) { 00144 // temp[dd] += b0[ii] * co[coefind*kdim + dd]; 00145 // } 00146 // coefind += 1; 00147 // } 00148 // for (int dd = 0; dd < kdim; ++dd) 00149 // restemp[dd] += temp[dd]*b1[jj]; 00150 // coefind += unum - uorder; 00151 // } 00152 00153 // // Copy from restemp to result 00154 // if (rational_) { 00155 // for (int dd = 0; dd < dim_; ++dd) { 00156 // result[dd] = restemp[dd]/restemp[kdim-1]; 00157 // } 00158 // } else { 00159 // for (int dd = 0; dd < dim_; ++dd) { 00160 // result[dd] = restemp[dd]; 00161 // } 00162 // } 00163 // } 00164 00165 00166 //=========================================================================== 00167 void 00168 SplineSurface::point(std::vector<Point>& result, double upar, double vpar, 00169 int derivs, bool u_from_right, bool v_from_right, 00170 double resolution) const 00171 //=========================================================================== 00172 { 00173 DEBUG_ERROR_IF(derivs < 0, "Negative number of derivatives makes no sense."); 00174 int totpts = (derivs + 1)*(derivs + 2)/2; 00175 int rsz = result.size(); 00176 DEBUG_ERROR_IF(rsz< totpts, "The vector of points must have sufficient size."); 00177 00178 for (int i = 0; i < totpts; ++i) { 00179 if (result[i].dimension() != dim_) { 00180 result[i].resize(dim_); 00181 } 00182 } 00183 00184 if (derivs == 0) { 00185 point(result[0], upar, vpar); 00186 return; 00187 } 00188 00189 // Take care of the rational case 00190 const std::vector<double>& co = rational_ ? rcoefs_ : coefs_; 00191 int kdim = dim_ + (rational_ ? 1 : 0); 00192 00193 // Make temporary storage for the basis values and a temporary 00194 // computation cache. 00195 Go::ScratchVect<double, 30> b0(basis_u_.order() * (derivs+1)); 00196 Go::ScratchVect<double, 30> b1(basis_v_.order() * (derivs+1)); 00197 Go::ScratchVect<double, 30> temp(kdim * totpts); 00198 Go::ScratchVect<double, 30> restemp(kdim * totpts); 00199 std::fill(restemp.begin(), restemp.end(), 0.0); 00200 // Compute the basis values and get some data about the spline spaces 00201 if (u_from_right) { 00202 basis_u_.computeBasisValues(upar, &b0[0], derivs, resolution); 00203 } else { 00204 basis_u_.computeBasisValuesLeft(upar, &b0[0], derivs, resolution); 00205 } 00206 int uleft = basis_u_.lastKnotInterval(); 00207 int uorder = basis_u_.order(); 00208 int unum = basis_u_.numCoefs(); 00209 if (v_from_right) { 00210 basis_v_.computeBasisValues(vpar, &b1[0], derivs, resolution); 00211 } else { 00212 basis_v_.computeBasisValuesLeft(vpar, &b1[0], derivs, resolution); 00213 } 00214 int vleft = basis_v_.lastKnotInterval(); 00215 int vorder = basis_v_.order(); 00216 // Compute the tensor product value 00217 int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); 00218 int derivs_plus1=derivs+1; 00219 for (int jj = 0; jj < vorder; ++jj) { 00220 int jjd=jj*(derivs_plus1); 00221 std::fill(temp.begin(), temp.end(), 0.0); 00222 00223 for (int ii = 0; ii < uorder; ++ii) { 00224 int iid=ii*(derivs_plus1); 00225 const double *co_p=&co[coefind*kdim]; 00226 for (int dd = 0; dd < kdim; ++dd,++co_p) { 00227 int temp_ind=dd; 00228 for (int vder = 0; vder < derivs_plus1; ++vder) { 00229 for (int uder = 0; uder < vder+1; ++uder) { 00230 temp[temp_ind] 00231 += b0[iid+vder - uder]*(*co_p); 00232 temp_ind+=kdim; 00233 } 00234 } 00235 } 00236 coefind += 1; 00237 } 00238 00239 for (int dd = 0; dd < kdim; ++dd) { 00240 int dercount = 0; 00241 for (int vder = 0; vder < derivs_plus1; ++vder) { 00242 for (int uder = 0; uder < vder + 1; ++uder) { 00243 restemp[dercount*kdim + dd] 00244 += temp[dercount*kdim + dd]*b1[uder + jjd]; 00245 ++dercount; 00246 } 00247 } 00248 } 00249 00250 coefind += unum - uorder; 00251 } 00252 // Copy from restemp to result 00253 if (rational_) { 00254 std::vector<double> restemp2(totpts*dim_); 00255 surface_ratder(&restemp[0], dim_, derivs, &restemp2[0]); 00256 for (int i = 0; i < totpts; ++i) { 00257 for (int dd = 0; dd < dim_; ++dd) { 00258 result[i][dd] = restemp2[i*dim_ + dd]; 00259 } 00260 } 00261 } else { 00262 double* restemp_it=restemp.begin(); 00263 for (int i = 0; i < totpts; ++i) { 00264 for (int dd = 0; dd < dim_; ++dd) { 00265 result[i][dd] = *restemp_it; 00266 ++restemp_it; 00267 } 00268 } 00269 } 00270 00271 } 00272 00273 #define NOT_FINISHED_YET 00274 #ifdef NOT_FINISHED_YET 00275 //=========================================================================== 00276 void SplineSurface::pointsGrid(int m1, int m2, int derivs, 00277 const double* basisvals1, 00278 const double* basisvals2, 00279 const int* knotint1, 00280 const int* knotint2, 00281 double* result, 00282 double* normals) const 00283 //=========================================================================== 00284 /* 00285 ********************************************************************* 00286 * 00287 ********************************************************************* 00288 * 00289 * PURPOSE : Evaluate this surface over an m1 * m2 00290 * grid, assuming that the B-splines have been 00291 * pre-evaluated, by s1504, over that grid. 00292 * The knots et1 and et2 and grid points (x[i],y[j]) are not needed. 00293 * Compute ider derivatives. 00294 * 00295 * INPUT : ps1 - Pointer to the surface to evaluate. 00296 * ider - Number of derivatives to calculate. 00297 * < 0 : No derivative calculated. 00298 * = 0 : Position calculated. 00299 * = 1 : Position and first derivative calculated. 00300 * etc. 00301 * m1 - Number of grid points in first direction. 00302 * m2 - Number of grid points in first direction. 00303 * ileft1 - Array of indexes to the intervals in the knotvector 00304 * in the first parameter direction where each subsequence 00305 * of k1 non-zero B-splines are located. 00306 * ileft2 - Array of indexes to the intervals in the knotvector 00307 * in the second parameter direction where each subsequence 00308 * of k2 non-zero B-splines are located. 00309 * ebder1 - Triple array of dimension [(ider+1)*k1*m1] containing 00310 * values of the k1 nonzero B-splines and their 00311 * derivatives at the points x[0],...,x[m1-1] 00312 * (i.e. pre-evaluated B-splines). 00313 * These numbers are stored in the following order: 00314 * First the (ider+1) derivatives of the first nonzero 00315 * B-spline at x[0]. Then the (ider+1) derivatives of 00316 * the second nonzero B-spline at x[0], etc. 00317 * Later we repeat for x[1],... etc. 00318 * 00319 * ebder2 - Triple array of dimension [(ider+1)*k2*m2] containing 00320 * values of the k2 nonzero B-splines and their 00321 * derivatives at the points y[0],...,y[m2-1] 00322 * 00323 * OUTPUT : eder - Array where the derivatives of the surface 00324 * are placed, dimension 00325 * idim * ((ider+1)(ider+2) / 2) * m1 * m2. 00326 * The sequence is position, 00327 * first derivative in first parameter direction, 00328 * first derivative in second parameter direction, 00329 * (2,0) derivative, (1,1) derivative, (0,2) 00330 * derivative, etc. at point (x[0],y[0]), 00331 * followed by the same information at (x[1],y[0]), 00332 * etc. 00333 * norm - Normals of surface. Is calculated if ider >= 1. 00334 * Dimension is idim*m1*m2. 00335 * The normals are normalized. @afr changed 00336 * jstat - status messages 00337 * = 2 : Surface is degenerate 00338 * at some point, normal 00339 * has zero length. 00340 * = 1 : Surface is close to 00341 * degenerate at some point 00342 * Angle between tangents, 00343 * less than angular tolerance. 00344 * = 0 : ok 00345 * < 0 : error 00346 * 00347 * METHOD : The code and method is similar to that of s1421 except that 00348 * the B-splines and their derivatives have already been 00349 * calculated (and are given in ebder1 and ebder2) and 00350 * we evaluate over a whole grid rather than at one point. 00351 * The method is to find the control points and control derivatives 00352 * of each isocurve in x (fixed y value). We then 00353 * evaluate at each x value along the isocurve. 00354 * 00355 * CALLS : s6err - errormacros.handling routine 00356 * s6strider - Make derivative of rational expression 00357 * 00358 * WRITTEN BY : Michael Floater, SINTEF, May 1998. 00359 ********************************************************************* 00360 */ 00361 { 00362 const double* ebder1 = basisvals1; 00363 const double* ebder2 = basisvals2; 00364 double* eder = result; 00365 double* norm = normals; 00366 const int* ileft1 = knotint1; 00367 const int* ileft2 = knotint2; 00368 int ider = derivs; 00369 int i1,i2; /* Loop variables. */ 00370 int i1pos,i2pos; /* Offset indexes. */ 00371 int kn1,kn2; /* The number of B-splines accociated with the knot 00372 vectors st1 and st2. */ 00373 int kk1,kk2; /* The polynomial order of the surface in the two 00374 directions. */ 00375 int kdim; /* The dimension of the space in which the surface 00376 lies. Equivalently, the number of components 00377 of each B-spline coefficient. */ 00378 int kleft1,kleft2; /* Local versions of knot intervals. */ 00379 int ki,kx,kjh; /* Control variables in for loops and for stepping 00380 through arrays. */ 00381 int kih2; /* Index for stepping through ebder2. */ 00382 int ky,kl,kl1,kl2; /* Control variables in for loops and for stepping 00383 through arrays. */ 00384 const double* scoef; /* The B-spline coefficients of the surface. 00385 This is an array of dimension [kn2*kn1*kdim]. */ 00386 double tt; /* Dummy variable used for holding an array element 00387 in a for loop. */ 00388 int size; /* Space occupied by points and derivs at one eval. */ 00389 int sizeh; /* Space occupied by homogeneous points and derivs . */ 00390 int size1,size2; /* Useful variables. */ 00391 /*int ederpos; */ /* Index of position in eder. */ 00392 /*int normpos; */ /* Index of position in norm. */ 00393 00394 int knumb2; /* Necessary size of ew */ 00395 00396 int tot,temp; /* Temporary variables. */ 00397 00398 00399 00400 kn1 = basis_u_.numCoefs(); 00401 kn2 = basis_v_.numCoefs(); 00402 kk1 = basis_u_.order(); 00403 kk2 = basis_v_.order(); 00404 00405 kdim = dim_; 00406 00407 if (rational_) { 00408 scoef = &rcoefs_[0]; 00409 kdim +=1; 00410 } else { 00411 scoef = &coefs_[0]; 00412 } 00413 00414 // @ afr: The following should have been declared outside this func. 00415 // Dvec eder(dimension() * ((ider+1)*(ider+2) / 2) * m1 * m2); 00416 // Dvec norm(dimension() * m1 * m2); 00417 00418 sizeh = kdim*(ider+1)*(ider+2)/2; 00419 double* sder = new double[sizeh]; 00420 double* enorm = new double[dim_]; 00421 00422 size = dim_*(ider+1)*(ider+2)/2; 00423 // int offset = m1*m2*kdim; 00424 size1 = (ider+1)*kk1; 00425 size2 = (ider+1)*kk2; 00426 00427 // Allocate space for B-spline values and derivatives and one work array. 00428 00429 knumb2 = kn1*(ider+1)*kdim; 00430 double* ew = new double[knumb2]; 00431 // double* ew_iterator = ew; // Iterator used for clearing ew. 00432 00433 double* norm_iterator = norm; 00434 double* eder_iterator = eder; 00435 00436 // Run through grid points in the y direction. 00437 for(i2=0, i2pos=0; i2<m2; i2++, i2pos += size2) { 00438 00439 kleft2 = ileft2[i2]; 00440 00441 /* Compute the control points (and control derivatives 00442 of the v = x[i2] isocurve. */ 00443 00444 /* Set all the elements of ew to 0. */ 00445 double* ew_iterator; 00446 for (ew_iterator=ew; ew_iterator<ew+knumb2; ew_iterator++) 00447 *ew_iterator = 0; 00448 00449 /* ki steps through the appropriate kk2 rows of B-spline coefficients 00450 while kih2 steps through the B-spline value and derivatives for the 00451 B-spline given by ki. */ 00452 00453 kih2 = i2pos; 00454 00455 for (ki=kleft2-kk2+1; ki<=kleft2; ki++) { 00456 /* kx counts through the ider+1 derivatives to be computed. 00457 kjh steps through ew once for each ki to accumulate the 00458 contribution from the different B-splines. 00459 kl1 points to the first component of the first B-spline 00460 coefficient in row no. ki of the B-spline coefficient matrix. 00461 */ 00462 00463 kjh = 0; kl1 = kdim*ki*kn1; 00464 for (kx=0; kx<=ider; kx++) { 00465 /* The value of the B-spline derivative is stored in tt while 00466 kl2 steps through the kdim components of all the B-spline 00467 coefficients that multiply nonzero B-splines along st1. 00468 */ 00469 00470 tt = ebder2[kih2++]; kl2 = kl1; 00471 for (kl=0; kl<kdim*kn1; kl++,kjh++,kl2++) { 00472 ew[kjh] += scoef[kl2]*tt; 00473 } 00474 } 00475 } 00476 00477 /* Run through grid points in the x direction 00478 evaluating along the iso-curve defined by y = y[i2]. */ 00479 for(i1=0, i1pos=0; i1<m1; i1++, i1pos += size1) { 00480 kleft1 = ileft1[i1]; 00481 00482 /* Set all the elements of sder to 0. */ 00483 for(ki=0; ki<sizeh; ki++) sder[ki] = 0; 00484 00485 for(ky=0; ky<=ider; ky++) { 00486 kl1 = kdim * (ky * kn1 + kleft1 - kk1 + 1); 00487 for(kx=0; kx<=ider-ky; kx++) { 00488 tot = kx + ky; 00489 temp = ((tot * (tot+1)) >> 1) + ky; 00490 kjh = temp * kdim; 00491 00492 for(ki=0; ki<kk1; ki++) { 00493 tt = ebder1[i1pos + (ider+1) * ki + kx]; 00494 for(kl=0; kl<kdim; kl++) { 00495 sder[kjh+kl] += ew[kl1 + kdim * ki + kl] * tt; 00496 } 00497 } 00498 } 00499 } 00500 00501 /* If rational surface calculate the derivatives based on 00502 derivatives in homogenous coordinates */ 00503 00504 // if (kind_ == 2 || kind_ == 4) { 00505 // // @ this must be implemented! 00506 // // s6strider(sder,ps1->idim,ider,eder+ederpos,&kstat); 00507 // } else { 00508 // for (ki=0; ki<sizeh; ki++) { 00509 // *(eder_iterator+ki) = sder[ki]; 00510 // } 00511 // } 00512 if (rational_) { 00513 surface_ratder(sder,dim_, ider, eder_iterator); 00514 // s6strider(sder,dim_,eder_iterator,&kstat); 00515 } else { 00516 for (ki=0; ki<sizeh; ki++) { 00517 *(eder_iterator+ki) = sder[ki]; 00518 } 00519 } 00520 00521 /* Calculate normal if idim==3 and ider>0. */ 00522 00523 if (ider>0 && kdim ==3) { 00524 // enorm = GoCrossProduct(eder_iterator + 3, 00525 // eder_iterator + 6); 00526 // GoNormalize(enorm); 00527 // copy(enorm,enorm+3,norm_iterator); 00528 double* i1 = eder + 3; 00529 double* i2 = eder + 6; 00530 norm_iterator[0] = i1[1]*i2[2] - i1[2]*i2[1]; 00531 norm_iterator[1] = i1[2]*i2[0] - i1[0]*i2[2]; 00532 norm_iterator[2] = i1[0]*i2[1] - i1[1]*i2[0]; 00533 double ssum = norm_iterator[0]*norm_iterator[0] 00534 + norm_iterator[1]*norm_iterator[1] 00535 + norm_iterator[2]*norm_iterator[2]; 00536 double invl = 1/sqrt(ssum); 00537 norm_iterator[0] = norm_iterator[0]*invl; 00538 norm_iterator[1] = norm_iterator[1]*invl; 00539 norm_iterator[2] = norm_iterator[2]*invl; 00540 } 00541 00542 eder_iterator += size; 00543 norm_iterator += kdim; 00544 } 00545 00546 } 00547 00548 00549 /* Free memory. */ 00550 00551 delete [] ew; 00552 delete [] sder; 00553 delete [] enorm; 00554 00555 00556 00557 /* Not enough memory. */ 00558 00559 /* kdim less than 1. */ 00560 00561 /* Polynomial order less than 1. */ 00562 00563 /* Fewer B-splines than the order. */ 00564 00565 /* Illegal derivative requested. */ 00566 00567 /* Error in lower level routine. */ 00568 00569 return; 00570 } 00571 #endif 00572 00573 00574 void SplineSurface::gridEvaluator(int num_u, int num_v, 00575 std::vector<double>& points, 00576 std::vector<double>& normals, 00577 std::vector<double>& param_u, 00578 std::vector<double>& param_v, 00579 bool normalize) const 00580 { 00581 ASSERT(dimension() == 3); 00582 ASSERT(num_u > 1 && num_v > 1); 00583 int kdim = rational_ ? dim_+1 : dim_; 00584 int ncoef_u = numCoefs_u(); 00585 // int ncoef_v = numCoefs_v(); 00586 int ord_u = order_u(); 00587 int ord_v = order_v(); 00588 00589 const double start_u = startparam_u(); 00590 const double start_v = startparam_v(); 00591 const double diff_u = endparam_u() - start_u; 00592 const double diff_v = endparam_v() - start_v; 00593 const double du=diff_u/(num_u-1); 00594 const double dv=diff_v/(num_v-1); 00595 00596 param_u.resize(num_u); 00597 param_v.resize(num_v); 00598 00599 for(int step_u = 0; step_u < num_u; ++step_u) 00600 param_u[step_u] = start_u+double(step_u)*du; 00601 for(int step_v = 0; step_v < num_v; ++step_v) 00602 param_v[step_v] = start_v+double(step_v)*dv; 00603 00604 vector<BasisDerivsSf> basis_func; 00605 00606 computeBasisGrid(param_u, param_v, basis_func); 00607 vector<double> der_u, der_v; 00608 points.resize(0); 00609 normals.resize(0); 00610 00611 vector<double> cf = rational_ ? rcoefs_ : coefs_; 00612 00613 double inv_weight = 1.0; 00614 vector<BasisDerivsSf>::const_iterator it = basis_func.begin(); 00615 for (int j = 0; j < num_v; ++j) 00616 { 00617 int coef_pos_v = ((*it).left_idx[1]-ord_v+1) * ncoef_u * kdim; 00618 for (int i = 0; i < num_u; ++i, ++it) 00619 { 00620 int coef_pos = coef_pos_v + ((*it).left_idx[0]-ord_u+1) * kdim; 00621 Point pt_coord(0.0, 0.0, 0.0); 00622 Point du_coord(0.0, 0.0, 0.0); 00623 Point dv_coord(0.0, 0.0, 0.0); 00624 int bas_pos = 0; 00625 // double denom = 0.0; 00626 for (int bas_v = 0; bas_v < ord_v; ++bas_v) 00627 { 00628 for (int bas_u = 0; bas_u < ord_u; ++bas_u, ++bas_pos) 00629 { 00630 if (rational_) 00631 inv_weight = 1.0 / cf[coef_pos + dim_]; 00632 for (int k = 0; k < dim_; ++k, ++coef_pos) 00633 { 00634 double current_coef = cf[coef_pos]; 00635 pt_coord[k] += current_coef * (*it).basisValues[bas_pos] * inv_weight; 00636 du_coord[k] += current_coef * (*it).basisDerivs_u[bas_pos] * inv_weight; 00637 dv_coord[k] += current_coef * (*it).basisDerivs_v[bas_pos] * inv_weight; 00638 } 00639 if (rational_) 00640 ++coef_pos; 00641 } 00642 coef_pos += (ncoef_u - ord_u)*kdim; 00643 } 00644 Point norm = du_coord % dv_coord; 00645 if (normalize) 00646 norm.normalize(); 00647 for (int k = 0; k < 3; ++k) 00648 { 00649 points.push_back(pt_coord[k]); 00650 normals.push_back(norm[k]); 00651 } 00652 } 00653 } 00654 } 00655 00656 // Same as above, but no normals. 00657 void SplineSurface::gridEvaluator(int num_u, int num_v, 00658 std::vector<double>& points, 00659 std::vector<double>& param_u, 00660 std::vector<double>& param_v, 00661 double start_u, double end_u, 00662 double start_v, double end_v) const 00663 { 00664 // ASSERT(dimension() == 3); 00665 ASSERT(num_u > 1 && num_v > 1); 00666 const double diff_u = end_u - start_u; 00667 const double diff_v = end_v - start_v; 00668 const double du=diff_u/(num_u-1); 00669 const double dv=diff_v/(num_v-1); 00670 00671 param_u.resize(num_u); 00672 param_v.resize(num_v); 00673 00674 for(int step_u = 0; step_u < num_u; ++step_u) 00675 param_u[step_u] = start_u+double(step_u)*du; 00676 for(int step_v = 0; step_v < num_v; ++step_v) 00677 param_v[step_v] = start_v+double(step_v)*dv; 00678 00679 gridEvaluator(points,param_u,param_v); 00680 } 00681 00682 // Added as separate method by kmo for usage in ICADA project. 00683 void SplineSurface::gridEvaluator(std::vector<double>& points, 00684 const std::vector<double>& param_u, 00685 const std::vector<double>& param_v) const 00686 { 00687 int num_u = param_u.size(); 00688 int num_v = param_v.size(); 00689 vector<double> basisvals_u(num_u * basis_u_.order()); 00690 vector<double> basisvals_v(num_v * basis_v_.order()); 00691 vector<int> knotinter_u(num_u * basis_u_.order()); 00692 vector<int> knotinter_v(num_v * basis_v_.order()); 00693 00694 basis_u_.computeBasisValues(¶m_u[0], ¶m_u[0]+param_u.size(), 00695 &basisvals_u[0], &knotinter_u[0], 0); 00696 basis_v_.computeBasisValues(¶m_v[0], ¶m_v[0]+param_v.size(), 00697 &basisvals_v[0], &knotinter_v[0], 0); 00698 00699 points.resize(num_u * num_v * dim_); 00700 pointsGrid(num_u, num_v, 0, 00701 &basisvals_u[0], &basisvals_v[0], 00702 &knotinter_u[0], &knotinter_v[0], 00703 &points[0], 0); 00704 } 00705 00706 // This thingie is intended to copy and adjust the existing 00707 // pointsGrid() method, so it would return correctly formatted results. 00708 void SplineSurface::pointsGridNoDerivs(int m1, int m2, 00709 const double* basisvals1, 00710 const double* basisvals2, 00711 const int* knotint1, 00712 const int* knotint2, 00713 double* result, 00714 double* normals, 00715 bool normalize) const 00716 { 00717 int i1,i2; /* Loop variables. */ 00718 int i1pos,i2pos; /* Offset indexes. */ 00719 int kn1,kn2; /* The number of B-splines accociated with the knot 00720 vectors st1 and st2. */ 00721 int kk1,kk2; /* The polynomial order of the surface in the two 00722 directions. */ 00723 int kleft1,kleft2; /* Local versions of knot intervals. */ 00724 int ki,kjh; /* Control variables in for loops and for stepping 00725 through arrays. */ 00726 int kih2; /* Index for stepping through basisvals2. */ 00727 int kl,kl1,kl2; /* Control variables in for loops and for stepping 00728 through arrays. */ 00729 const double* scoef; /* The B-spline coefficients of the surface. 00730 This is an array of dimension [kn2*kn1*kdim]. */ 00731 double tt; /* Dummy variable used for holding an array element 00732 in a for loop. */ 00733 int size1,size2; /* Useful variables. */ 00734 00735 int knumb2; /* Necessary size of ew */ 00736 00737 00738 00739 00740 kn1 = basis_u_.numCoefs(); 00741 kn2 = basis_v_.numCoefs(); 00742 kk1 = basis_u_.order(); 00743 kk2 = basis_v_.order(); 00744 00745 // Check that dim_=3 && rational_=false; 00746 00747 scoef = &coefs_[0]; 00748 00749 00750 ScratchVect<double, 9> sder(9); 00751 ScratchVect<double, 3> enorm(3); 00752 00753 // int offset = m1*m2*9; 00754 size1 = 2*kk1; 00755 size2 = 2*kk2; 00756 00757 // Allocate space for B-spline values and derivatives and one work array. 00758 00759 knumb2 = kn1*6; 00760 ScratchVect<double, 300> ew(knumb2); 00761 00762 double* normals_iterator = normals; 00763 double* result_iterator = result; 00764 00765 // Run through grid points in the y direction. 00766 for(i2=0, i2pos=0; i2<m2; i2++, i2pos += size2) { 00767 00768 kleft2 = knotint2[i2]; 00769 00770 /* Compute the control points (and control derivatives 00771 of the v = x[i2] isocurve. */ 00772 00773 /* Set all the elements of ew to 0. */ 00774 for (int i = 0; i<knumb2; i++) 00775 ew[i] = 0; 00776 00777 /* ki steps through the appropriate kk2 rows of B-spline coefficients 00778 while kih2 steps through the B-spline value and derivatives for the 00779 B-spline given by ki. */ 00780 00781 kih2 = i2pos; 00782 00783 for (ki=kleft2-kk2+1; ki<=kleft2; ki++) { 00784 /* kx counts through the ider+1 derivatives to be computed. 00785 kjh steps through ew once for each ki to accumulate the 00786 contribution from the different B-splines. 00787 kl1 points to the first component of the first B-spline 00788 coefficient in row no. ki of the B-spline coefficient matrix. 00789 */ 00790 00791 kjh = 0; kl1 = 3*ki*kn1; 00792 /* The value of the B-spline derivative is stored in tt while 00793 kl2 steps through the kdim components of all the B-spline 00794 coefficients that multiply nonzero B-splines along st1. 00795 */ 00796 tt = basisvals2[kih2++]; kl2 = kl1; 00797 for (kl=0; kl<3*kn1; kl++,kjh++,kl2++) { 00798 ew[kjh] += scoef[kl2]*tt; 00799 } 00800 tt = basisvals2[kih2++]; kl2 = kl1; 00801 for (kl=0; kl<3*kn1; kl++,kjh++,kl2++) { 00802 ew[kjh] += scoef[kl2]*tt; 00803 } 00804 } 00805 00806 /* Run through grid points in the x direction 00807 evaluating along the iso-curve defined by y = y[i2]. */ 00808 for(i1=0, i1pos=0; i1<m1; i1++, i1pos += size1) { 00809 kleft1 = knotint1[i1]; 00810 00811 /* Set all the elements of sder to 0. */ 00812 sder[0] = sder[1] = sder[2] = sder[3] = sder[4] = 0; 00813 sder[5] = sder[6] = sder[7] = sder[8] = 0; 00814 00815 kl1 = 3 * (kleft1 - kk1 + 1); 00816 for(ki=0; ki<kk1; ki++) { 00817 tt = basisvals1[i1pos + ki*2]; 00818 sder[0] += ew[kl1 + ki*3 + 0] * tt; 00819 sder[1] += ew[kl1 + ki*3 + 1] * tt; 00820 sder[2] += ew[kl1 + ki*3 + 2] * tt; 00821 } 00822 for(ki=0; ki<kk1; ki++) { 00823 tt = basisvals1[i1pos + ki*2 + 1]; 00824 sder[3] += ew[kl1 + ki*3 + 0] * tt; 00825 sder[4] += ew[kl1 + ki*3 + 1] * tt; 00826 sder[5] += ew[kl1 + ki*3 + 2] * tt; 00827 } 00828 00829 kl1 = 3 * (kn1 + kleft1 - kk1 + 1); 00830 for(ki=0; ki<kk1; ki++) { 00831 tt = basisvals1[i1pos + ki*2]; 00832 sder[6] += ew[kl1 + ki*3 + 0] * tt; 00833 sder[7] += ew[kl1 + ki*3 + 1] * tt; 00834 sder[8] += ew[kl1 + ki*3 + 2] * tt; 00835 } 00836 kjh = 6; 00837 00838 *(result_iterator++) = sder[0]; 00839 *(result_iterator++) = sder[1]; 00840 *(result_iterator++) = sder[2]; 00841 00842 double* i1 = &sder[3]; 00843 double* i2 = &sder[6]; 00844 normals_iterator[0] = i1[1]*i2[2] - i1[2]*i2[1]; 00845 normals_iterator[1] = i1[2]*i2[0] - i1[0]*i2[2]; 00846 normals_iterator[2] = i1[0]*i2[1] - i1[1]*i2[0]; 00847 if (normalize) 00848 { 00849 double ssum = normals_iterator[0]*normals_iterator[0] 00850 + normals_iterator[1]*normals_iterator[1] 00851 + normals_iterator[2]*normals_iterator[2]; 00852 double invl = 1/sqrt(ssum); 00853 normals_iterator[0] = normals_iterator[0]*invl; 00854 normals_iterator[1] = normals_iterator[1]*invl; 00855 normals_iterator[2] = normals_iterator[2]*invl; 00856 } 00857 normals_iterator += 3; 00858 } 00859 00860 } 00861 00862 return; 00863 } 00864 00865 00866 void SplineSurface::computeBasis(double param[], 00867 std::vector< double > &basisValues, 00868 std::vector< double > &basisDerivs_u, 00869 std::vector< double > &basisDerivs_v) const 00870 { 00871 int kk1 = basis_u_.order(); 00872 int kk2 = basis_v_.order(); 00873 int nn1 = basis_u_.numCoefs(); 00874 vector<double> basisvals_u(2*kk1); 00875 vector<double> basisvals_v(2*kk2); 00876 00877 // Compute basis values 00878 basis_u_.computeBasisValues(param[0], &basisvals_u[0], 1); 00879 basis_v_.computeBasisValues(param[1], &basisvals_v[0], 1); 00880 00881 // Accumulate 00882 int ki, kj, kr; 00883 basisValues.resize(kk1*kk2); 00884 basisDerivs_u.resize(kk1*kk2); 00885 basisDerivs_v.resize(kk1*kk2); 00886 vector<double> weights(kk1*kk2); 00887 if (rational_) 00888 { 00889 int kdim = dim_ + 1; 00890 int uleft = basis_u_.lastKnotInterval() - kk1 + 1; 00891 int vleft = basis_v_.lastKnotInterval() - kk2 + 1; 00892 for (kj=vleft, kr=0; kj<vleft+kk2; ++kj) 00893 for (ki=uleft; ki<uleft+kk1; ++ki) 00894 weights[kr++] = rcoefs_[(kj*nn1+ki)*kdim+dim_]; 00895 } 00896 00897 accumulateBasis(&basisvals_u[0], kk1, &basisvals_v[0], kk2, 00898 rational_ ? &weights[0] : NULL, &basisValues[0], 00899 &basisDerivs_u[0], &basisDerivs_v[0]); 00900 00901 } 00902 00903 //=========================================================================== 00904 void SplineSurface::computeBasisGrid(const Dvector& param_u, 00905 const Dvector& param_v, 00906 Dmatrix& basisValues) const 00907 //=========================================================================== 00908 { 00909 int numu = (int)param_u.size(); 00910 int numv = (int)param_v.size(); 00911 int uorder = basis_u_.order(); 00912 int vorder = basis_v_.order(); 00913 int ucoefs = basis_u_.numCoefs(); 00914 int vcoefs = basis_v_.numCoefs(); 00915 00916 vector<double> basisvals_u(numu * uorder); 00917 vector<double> basisvals_v(numv * vorder); 00918 vector<int> left_u(numu); 00919 vector<int> left_v(numv); 00920 00921 // Compute basis values 00922 basis_u_.computeBasisValues(¶m_u[0], ¶m_u[0]+param_u.size(), 00923 &basisvals_u[0], &left_u[0]); 00924 basis_v_.computeBasisValues(¶m_v[0], ¶m_v[0]+param_v.size(), 00925 &basisvals_v[0], &left_v[0]); 00926 00927 // Initiate to zero 00928 int ki; 00929 int num_coefs = ucoefs*vcoefs; 00930 int num_par = numu*numv; 00931 basisValues.resize(num_par); 00932 for (ki=0; ki<num_par; ++ki) 00933 basisValues[ki].assign(num_coefs, 0.0); 00934 00935 // Fetch all weights 00936 vector<double> weights, currw; 00937 if (rational_) 00938 { 00939 currw.resize(uorder*vorder); 00940 weights.resize(ucoefs*vcoefs); 00941 getWeights(weights); 00942 } 00943 00944 // For all points 00945 int kj, kh, kv; 00946 int idx1, idx2, idx4, idx5; 00947 int numorder = uorder*vorder; 00948 vector<double> tmpVal(numorder); 00949 for (kh=0, kj=0, idx2=0; kj<numv; ++kj, idx2+=vorder) 00950 { 00951 for (ki=0, idx1=0; ki<numu; ++ki, ++kh, idx1+=uorder) 00952 { 00953 int uleft = left_u[ki] - uorder + 1; 00954 int vleft = left_v[kj] - vorder + 1; 00955 if (rational_) 00956 { 00957 // Collect relevant weights 00958 vector<double>::iterator wgt = weights.begin() + vleft*ucoefs; 00959 vector<double>::iterator currwgt = currw.begin(); 00960 for (kv=0; kv<vorder; ++kv, wgt+=ucoefs, currwgt+=uorder) 00961 { 00962 std::copy(wgt+uleft, wgt+uleft+uorder, currwgt); 00963 } 00964 } 00965 00966 accumulateBasis(&basisvals_u[idx1], uorder, &basisvals_v[idx2], 00967 vorder, rational_ ? &currw[0] : NULL, &tmpVal[0]); 00968 00969 // Copy results into output array 00970 for (kv=0, idx4=0, idx5=vleft*ucoefs+uleft; kv<vorder; ++kv, idx4+=uorder, idx5+=ucoefs) 00971 { 00972 std::copy(tmpVal.begin()+idx4, tmpVal.begin()+idx4+uorder, 00973 basisValues[kh].begin()+idx5); 00974 } 00975 } 00976 } 00977 00978 } 00979 00980 //=========================================================================== 00981 void SplineSurface::computeBasisGrid(const Dvector& param_u, 00982 const Dvector& param_v, 00983 Dmatrix& basisValues, 00984 Dmatrix& basisDerivs_u, 00985 Dmatrix& basisDerivs_v) const 00986 //=========================================================================== 00987 { 00988 int derivs = 1; // Compute position and 1. derivative 00989 int numu = (int)param_u.size(); 00990 int numv = (int)param_v.size(); 00991 int uorder = basis_u_.order(); 00992 int vorder = basis_v_.order(); 00993 int ucoefs = basis_u_.numCoefs(); 00994 int vcoefs = basis_v_.numCoefs(); 00995 00996 vector<double> basisvals_u(numu * uorder * (derivs + 1)); 00997 vector<double> basisvals_v(numv * vorder * (derivs + 1)); 00998 vector<int> left_u(numu); 00999 vector<int> left_v(numv); 01000 01001 // Compute basis values 01002 basis_u_.computeBasisValues(¶m_u[0], ¶m_u[0]+param_u.size(), 01003 &basisvals_u[0], &left_u[0], derivs); 01004 basis_v_.computeBasisValues(¶m_v[0], ¶m_v[0]+param_v.size(), 01005 &basisvals_v[0], &left_v[0], derivs); 01006 01007 // Initiate to zero 01008 int ki; 01009 int num_coefs = ucoefs*vcoefs; 01010 int num_par = numu*numv; 01011 basisValues.resize(num_par); 01012 basisDerivs_u.resize(num_par); 01013 basisDerivs_v.resize(num_par); 01014 for (ki=0; ki<num_par; ++ki) 01015 { 01016 basisValues[ki].assign(num_coefs, 0.0); 01017 basisDerivs_u[ki].assign(num_coefs, 0.0); 01018 basisDerivs_v[ki].assign(num_coefs, 0.0); 01019 } 01020 01021 // Fetch all weights 01022 vector<double> weights, currw; 01023 if (rational_) 01024 { 01025 currw.resize(uorder*vorder); 01026 weights.resize(ucoefs*vcoefs); 01027 getWeights(weights); 01028 } 01029 01030 // For all points 01031 int kj, kh, kv; 01032 int idx1, idx2, idx4, idx5; 01033 int numorder = uorder*vorder; 01034 vector<double> tmpVal(numorder), tmpDer_u(numorder), tmpDer_v(numorder); 01035 for (kh=0, kj=0, idx2=0; kj<numv; ++kj, idx2+=2*vorder) 01036 { 01037 for (ki=0, idx1=0; ki<numu; ++ki, ++kh, idx1+=2*uorder) 01038 { 01039 int uleft = left_u[ki] - uorder + 1; 01040 int vleft = left_v[kj] - vorder + 1; 01041 if (rational_) 01042 { 01043 // Collect relevant weights 01044 vector<double>::iterator wgt = weights.begin() + vleft*ucoefs; 01045 vector<double>::iterator currwgt = currw.begin(); 01046 for (kv=0; kv<vorder; ++kv, wgt+=ucoefs, currwgt+=uorder) 01047 { 01048 std::copy(wgt+uleft, wgt+uleft+uorder, currwgt); 01049 } 01050 } 01051 01052 accumulateBasis(&basisvals_u[idx1], uorder, &basisvals_v[idx2], 01053 vorder, rational_ ? &currw[0] : NULL, 01054 &tmpVal[0], &tmpDer_u[0], &tmpDer_v[0]); 01055 01056 // Copy results into output array 01057 for (kv=0, idx4=0, idx5=vleft*ucoefs+uleft; kv<vorder; ++kv, idx4+=uorder, idx5+=ucoefs) 01058 { 01059 std::copy(tmpVal.begin()+idx4, tmpVal.begin()+idx4+uorder, 01060 basisValues[kh].begin()+idx5); 01061 std::copy(tmpDer_u.begin()+idx4, tmpDer_u.begin()+idx4+uorder, 01062 basisDerivs_u[kh].begin()+idx5); 01063 std::copy(tmpDer_v.begin()+idx4, tmpDer_v.begin()+idx4+uorder, 01064 basisDerivs_v[kh].begin()+idx5); 01065 } 01066 } 01067 } 01068 01069 } 01070 01071 //=========================================================================== 01072 void SplineSurface::computeBasisGrid(const Dvector& param_u, 01073 const Dvector& param_v, 01074 vector<BasisPtsSf>& result) const 01075 //=========================================================================== 01076 { 01077 int numu = (int)param_u.size(); 01078 int numv = (int)param_v.size(); 01079 int uorder = basis_u_.order(); 01080 int vorder = basis_v_.order(); 01081 int ucoefs = basis_u_.numCoefs(); 01082 int vcoefs = basis_v_.numCoefs(); 01083 01084 vector<double> basisvals_u(numu * uorder); 01085 vector<double> basisvals_v(numv * vorder); 01086 vector<int> left_u(numu); 01087 vector<int> left_v(numv); 01088 01089 // Compute basis values 01090 basis_u_.computeBasisValues(¶m_u[0], ¶m_u[0]+param_u.size(), 01091 &basisvals_u[0], &left_u[0]); 01092 basis_v_.computeBasisValues(¶m_v[0], ¶m_v[0]+param_v.size(), 01093 &basisvals_v[0], &left_v[0]); 01094 01095 // Initiate output 01096 result.resize(numu*numv); 01097 01098 // Fetch all weights 01099 vector<double> weights, currw; 01100 if (rational_) 01101 { 01102 currw.resize(uorder*vorder); 01103 weights.resize(ucoefs*vcoefs); 01104 getWeights(weights); 01105 } 01106 01107 // For all points 01108 int ki, kj, kh, kv; 01109 int idx1, idx2; 01110 for (kh=0, kj=0, idx2=0; kj<numv; ++kj, idx2+=vorder) 01111 { 01112 for (ki=0, idx1=0; ki<numu; ++ki, ++kh, idx1+=uorder) 01113 { 01114 result[kh].preparePts(param_u[ki], param_v[kj], 01115 left_u[ki], left_v[kj], 01116 uorder*vorder); 01117 01118 01119 if (rational_) 01120 { 01121 // Collect relevant weights 01122 int uleft = left_u[ki] - uorder + 1; 01123 int vleft = left_v[kj] - vorder + 1; 01124 vector<double>::iterator wgt = weights.begin() + vleft*ucoefs; 01125 vector<double>::iterator currwgt = currw.begin(); 01126 for (kv=0; kv<vorder; ++kv, wgt+=ucoefs, currwgt+=uorder) 01127 { 01128 std::copy(wgt+uleft, wgt+uleft+uorder, currwgt); 01129 } 01130 } 01131 01132 accumulateBasis(&basisvals_u[idx1], uorder, &basisvals_v[idx2], 01133 vorder, rational_ ? &currw[0] : NULL, 01134 &result[kh].basisValues[0]); 01135 } 01136 } 01137 } 01138 01139 01140 //=========================================================================== 01141 void SplineSurface::computeBasisGrid(const Dvector& param_u, 01142 const Dvector& param_v, 01143 vector<BasisDerivsSf>& result) const 01144 //=========================================================================== 01145 { 01146 int derivs = 1; // Compute position and 1. derivative 01147 int numu = (int)param_u.size(); 01148 int numv = (int)param_v.size(); 01149 int uorder = basis_u_.order(); 01150 int vorder = basis_v_.order(); 01151 int ucoefs = basis_u_.numCoefs(); 01152 int vcoefs = basis_v_.numCoefs(); 01153 01154 vector<double> basisvals_u(numu * uorder * (derivs + 1)); 01155 vector<double> basisvals_v(numv * vorder * (derivs + 1)); 01156 vector<int> left_u(numu); 01157 vector<int> left_v(numv); 01158 01159 // Compute basis values 01160 basis_u_.computeBasisValues(¶m_u[0], ¶m_u[0]+param_u.size(), 01161 &basisvals_u[0], &left_u[0], derivs); 01162 basis_v_.computeBasisValues(¶m_v[0], ¶m_v[0]+param_v.size(), 01163 &basisvals_v[0], &left_v[0], derivs); 01164 01165 // Initiate output 01166 result.resize(numu*numv); 01167 01168 // Fetch all weights 01169 vector<double> weights, currw; 01170 if (rational_) 01171 { 01172 currw.resize(uorder*vorder); 01173 weights.resize(ucoefs*vcoefs); 01174 getWeights(weights); 01175 } 01176 01177 // For all points 01178 int ki, kj, kh, kv; 01179 int idx1, idx2; 01180 for (kh=0, kj=0, idx2=0; kj<numv; ++kj, idx2+=2*vorder) 01181 { 01182 for (ki=0, idx1=0; ki<numu; ++ki, ++kh, idx1+=2*uorder) 01183 { 01184 result[kh].prepareDerivs(param_u[ki], param_v[kj], 01185 left_u[ki], left_v[kj], 01186 uorder*vorder); 01187 01188 01189 if (rational_) 01190 { 01191 // Collect relevant weights 01192 int uleft = left_u[ki] - uorder + 1; 01193 int vleft = left_v[kj] - vorder + 1; 01194 vector<double>::iterator wgt = weights.begin() + vleft*ucoefs; 01195 vector<double>::iterator currwgt = currw.begin(); 01196 for (kv=0; kv<vorder; ++kv, wgt+=ucoefs, currwgt+=uorder) 01197 { 01198 std::copy(wgt+uleft, wgt+uleft+uorder, currwgt); 01199 } 01200 } 01201 01202 accumulateBasis(&basisvals_u[idx1], uorder, &basisvals_v[idx2], 01203 vorder, rational_ ? &currw[0] : NULL, 01204 &result[kh].basisValues[0], 01205 &result[kh].basisDerivs_u[0], &result[kh].basisDerivs_v[0]); 01206 } 01207 } 01208 } 01209 01210 01211 //=========================================================================== 01212 void SplineSurface::computeBasisGrid(const Dvector& param_u, 01213 const Dvector& param_v, 01214 vector<BasisDerivsSf2>& result) const 01215 //=========================================================================== 01216 { 01217 int derivs = 2; // Compute position, 1. and 2. derivative 01218 int numu = (int)param_u.size(); 01219 int numv = (int)param_v.size(); 01220 int uorder = basis_u_.order(); 01221 int vorder = basis_v_.order(); 01222 int ucoefs = basis_u_.numCoefs(); 01223 int vcoefs = basis_v_.numCoefs(); 01224 01225 vector<double> basisvals_u(numu * uorder * (derivs + 1)); 01226 vector<double> basisvals_v(numv * vorder * (derivs + 1)); 01227 vector<int> left_u(numu); 01228 vector<int> left_v(numv); 01229 01230 // Compute basis values 01231 basis_u_.computeBasisValues(¶m_u[0], ¶m_u[0]+param_u.size(), 01232 &basisvals_u[0], &left_u[0], derivs); 01233 basis_v_.computeBasisValues(¶m_v[0], ¶m_v[0]+param_v.size(), 01234 &basisvals_v[0], &left_v[0], derivs); 01235 01236 // Initiate output 01237 result.resize(numu*numv); 01238 01239 // Fetch all weights 01240 vector<double> weights, currw; 01241 if (rational_) 01242 { 01243 currw.resize(uorder*vorder); 01244 weights.resize(ucoefs*vcoefs); 01245 getWeights(weights); 01246 } 01247 01248 // For all points 01249 int ki, kj, kh, kv; 01250 int idx1, idx2; 01251 for (kh=0, kj=0, idx2=0; kj<numv; ++kj, idx2+=3*vorder) 01252 { 01253 for (ki=0, idx1=0; ki<numu; ++ki, ++kh, idx1+=3*uorder) 01254 { 01255 result[kh].prepareDerivs(param_u[ki], param_v[kj], 01256 left_u[ki], left_v[kj], 01257 uorder*vorder); 01258 01259 if (rational_) 01260 { 01261 // Collect relevant weights 01262 int uleft = left_u[ki] - uorder + 1; 01263 int vleft = left_v[kj] - vorder + 1; 01264 vector<double>::iterator wgt = weights.begin() + vleft*ucoefs; 01265 vector<double>::iterator currwgt = currw.begin(); 01266 for (kv=0; kv<vorder; ++kv, wgt+=ucoefs, currwgt+=uorder) 01267 { 01268 std::copy(wgt+uleft, wgt+uleft+uorder, currwgt); 01269 } 01270 } 01271 01272 accumulateBasis(&basisvals_u[idx1], uorder, &basisvals_v[idx2], 01273 vorder, rational_ ? &currw[0] : NULL, 01274 &result[kh].basisValues[0], 01275 &result[kh].basisDerivs_u[0], &result[kh].basisDerivs_v[0], 01276 &result[kh].basisDerivs_uu[0], &result[kh].basisDerivs_uv[0], &result[kh].basisDerivs_vv[0]); 01277 } 01278 } 01279 01280 } 01281 01282 01283 01284 //=========================================================================== 01285 void SplineSurface::accumulateBasis(double* basisvals_u, int uorder, 01286 double* basisvals_v, int vorder, 01287 double* weights, 01288 double* basisValues) const 01289 //=========================================================================== 01290 { 01291 int ki, kj, kr; 01292 if (rational_) 01293 { 01294 // Compute denominator and derivatives thereof 01295 double sum = 0.0; 01296 for (kj=0, kr=0; kj<vorder; ++kj) 01297 for (ki=0; ki<uorder; ++ki, ++kr) 01298 sum += basisvals_u[ki]*basisvals_v[kj]*weights[kr]; 01299 01300 // Compute rational expression 01301 // double sum2 = sum*sum; 01302 for (kj=0, kr=0; kj<vorder; ++kj) 01303 for (ki=0; ki<uorder; ++ki, ++kr) 01304 basisValues[kr] = basisvals_u[ki]*basisvals_v[kj]*weights[kr]/sum; 01305 } 01306 else 01307 { 01308 // Multiply basis values in the three parameter directions 01309 for (kj=0, kr=0; kj<vorder; ++kj) 01310 for (ki=0; ki<uorder; ++ki, ++kr) 01311 basisValues[kr] = basisvals_u[ki]*basisvals_v[kj]; 01312 } 01313 01314 } 01315 01316 //=========================================================================== 01317 void SplineSurface::accumulateBasis(double* basisvals_u, int uorder, 01318 double* basisvals_v, int vorder, 01319 double* weights, 01320 double* basisValues, 01321 double* basisDerivs_u, 01322 double* basisDerivs_v) const 01323 //=========================================================================== 01324 { 01325 int ki, kj, kr; 01326 if (rational_) 01327 { 01328 // Compute denominator and derivatives thereof 01329 double sum = 0.0, dusum = 0.0, dvsum = 0.0; 01330 for (kj=0, kr=0; kj<vorder; ++kj) 01331 for (ki=0; ki<uorder; ++ki, ++kr) 01332 { 01333 sum += basisvals_u[2*ki]*basisvals_v[2*kj]*weights[kr]; 01334 dusum += basisvals_u[2*ki+1]*basisvals_v[2*kj]*weights[kr]; 01335 dvsum += basisvals_u[2*ki]*basisvals_v[2*kj+1]*weights[kr]; 01336 } 01337 01338 // Compute rational expression 01339 double sum2 = sum*sum; 01340 for (kj=0, kr=0; kj<vorder; ++kj) 01341 { 01342 double tmp2 = (basisvals_v[2*kj+1]*sum - basisvals_v[2*kj]*dvsum)/sum2; 01343 double fac = basisvals_v[2*kj]/sum2; 01344 for (ki=0; ki<uorder; ++ki, ++kr) 01345 { 01346 basisValues[kr] = basisvals_u[2*ki]*basisvals_v[2*kj]*weights[kr]/sum; 01347 double tmp1 = (basisvals_u[2*ki+1]*sum - basisvals_u[2*ki]*dusum)*weights[kr]*fac; 01348 basisDerivs_u[kr] = tmp1; 01349 basisDerivs_v[kr] = tmp2*weights[kr]*basisvals_u[2*ki]; 01350 } 01351 } 01352 } 01353 else 01354 { 01355 // Multiply basis values in the two parameter directions 01356 for (kj=0, kr=0; kj<vorder; ++kj) 01357 for (ki=0; ki<uorder; ++ki, ++kr) 01358 { 01359 basisValues[kr] = basisvals_u[2*ki]*basisvals_v[2*kj]; 01360 basisDerivs_u[kr] = basisvals_u[2*ki+1]*basisvals_v[2*kj]; 01361 basisDerivs_v[kr] = basisvals_u[2*ki]*basisvals_v[2*kj+1]; 01362 } 01363 } 01364 01365 } 01366 01367 01368 //=========================================================================== 01369 void SplineSurface::accumulateBasis(double* basisvals_u, int uorder, 01370 double* basisvals_v, int vorder, 01371 double* weights, 01372 double* basisValues, 01373 double* basisDerivs_u, 01374 double* basisDerivs_v, 01375 double* basisDerivs_uu, 01376 double* basisDerivs_uv, 01377 double* basisDerivs_vv) const 01378 //=========================================================================== 01379 { 01380 int ki, kj, kr; 01381 if (rational_) 01382 { 01383 // Compute denominator and derivatives thereof 01384 double 01385 sum = 0.0, dusum = 0.0, dvsum = 0.0, 01386 duusum = 0.0, duvsum = 0.0, dvvsum = 0.0; 01387 for (kj=0, kr=0; kj<vorder; ++kj) 01388 for (ki=0; ki<uorder; ++ki, ++kr) 01389 { 01390 sum += basisvals_u[3*ki]*basisvals_v[3*kj]*weights[kr]; 01391 dusum += basisvals_u[3*ki+1]*basisvals_v[3*kj]*weights[kr]; 01392 dvsum += basisvals_u[3*ki]*basisvals_v[3*kj+1]*weights[kr]; 01393 duusum += basisvals_u[3*ki+2]*basisvals_v[3*kj]*weights[kr]; 01394 duvsum += basisvals_u[3*ki+1]*basisvals_v[3*kj+1]*weights[kr]; 01395 dvvsum += basisvals_u[3*ki]*basisvals_v[3*kj+2]*weights[kr]; 01396 } 01397 01398 // Compute rational expression 01399 double sum2 = sum*sum; 01400 double sum3 = sum*sum2; 01401 for (kj=0, kr=0; kj<vorder; ++kj) 01402 { 01403 double tmp_v_v = (basisvals_v[3*kj+1]*sum - basisvals_v[3*kj]*dvsum)/sum2; 01404 double tmp_v_vv = (basisvals_v[3*kj+2]*sum2 + 2.0*basisvals_v[3*kj]*dvsum*dvsum 01405 - basisvals_v[3*kj]*dvvsum*sum - 2.0*basisvals_v[3*kj+1]*dvsum*sum) / sum3; 01406 double fac_v2 = basisvals_v[3*kj] / sum2; 01407 double fac_v3 = basisvals_v[3*kj] / sum3; 01408 for (ki=0; ki<uorder; ++ki, ++kr) 01409 { 01410 01411 basisValues[kr] = basisvals_u[3*ki]*basisvals_v[3*kj]*weights[kr]/sum; 01412 01413 basisDerivs_u[kr] = (basisvals_u[3*ki+1]*sum - basisvals_u[3*ki]*dusum) * weights[kr] * fac_v2; 01414 basisDerivs_v[kr] = tmp_v_v * weights[kr] * basisvals_u[3*ki]; 01415 01416 basisDerivs_uu[kr] = (basisvals_u[3*ki+2]*sum2 + 2.0*basisvals_u[3*ki]*dusum*dusum 01417 - basisvals_u[3*ki]*duusum*sum - 2.0*basisvals_u[3*ki+1]*dusum*sum) * weights[kr] * fac_v3; 01418 basisDerivs_uv[kr] = (basisvals_u[3*ki+1]*basisvals_v[3*kj+1]*sum2 + 2.0*basisvals_u[3*ki]*basisvals_v[3*kj]*dusum*dvsum 01419 - basisvals_u[3*ki+1]*basisvals_v[3*kj]*dvsum*sum - basisvals_u[3*ki]*basisvals_v[3*kj+1]*dusum*sum 01420 - basisvals_u[3*ki]*basisvals_v[3*kj]*duvsum*sum) * weights[kr] / sum3; 01421 basisDerivs_vv[kr] = tmp_v_vv * weights[kr] * basisvals_u[3*ki]; 01422 } 01423 } 01424 } 01425 else 01426 { 01427 // Multiply basis values in the two parameter directions 01428 for (kj=0, kr=0; kj<vorder; ++kj) 01429 for (ki=0; ki<uorder; ++ki, ++kr) 01430 { 01431 basisValues[kr] = basisvals_u[3*ki]*basisvals_v[3*kj]; 01432 basisDerivs_u[kr] = basisvals_u[3*ki+1]*basisvals_v[3*kj]; 01433 basisDerivs_v[kr] = basisvals_u[3*ki]*basisvals_v[3*kj+1]; 01434 basisDerivs_uu[kr] = basisvals_u[3*ki+2]*basisvals_v[3*kj]; 01435 basisDerivs_uv[kr] = basisvals_u[3*ki+1]*basisvals_v[3*kj+1]; 01436 basisDerivs_vv[kr] = basisvals_u[3*ki]*basisvals_v[3*kj+2]; 01437 } 01438 } 01439 01440 } 01441 01442 01443 } // namespace Go 01444 01445 01446