?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed: gcc8 base update
[opensg.git] / Source / System / NodeCores / Drawables / Nurbs / Internal / OSGBSplineTensorSurface.cpp
blob8609d89c87801bbb0b4c7db45e5ade01b144311d
1 /*---------------------------------------------------------------------------*\
2 * OpenSG NURBS Library *
3 * *
4 * *
5 * Copyright (C) 2001-2006 by the University of Bonn, Computer Graphics Group*
6 * *
7 * http://cg.cs.uni-bonn.de/ *
8 * *
9 * contact: edhellon@cs.uni-bonn.de, guthe@cs.uni-bonn.de, rk@cs.uni-bonn.de *
10 * *
11 \*---------------------------------------------------------------------------*/
12 /*---------------------------------------------------------------------------*\
13 * License *
14 * *
15 * This library is free software; you can redistribute it and/or modify it *
16 * under the terms of the GNU Library General Public License as published *
17 * by the Free Software Foundation, version 2. *
18 * *
19 * This library is distributed in the hope that it will be useful, but *
20 * WITHOUT ANY WARRANTY; without even the implied warranty of *
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
22 * Library General Public License for more details. *
23 * *
24 * You should have received a copy of the GNU Library General Public *
25 * License along with this library; if not, write to the Free Software *
26 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *
27 * *
28 \*---------------------------------------------------------------------------*/
29 /*---------------------------------------------------------------------------*\
30 * Changes *
31 * *
32 * *
33 * *
34 * *
35 * *
36 * *
37 \*---------------------------------------------------------------------------*/
38 #ifdef WIN32
39 # pragma warning (disable : 985)
40 #endif
41 #include "OSGBSplineTensorSurface.h"
42
43 OSG_USING_NAMESPACE
44
45
46
47 // FIXME: clean up all of the copy'n'paste mess...
48 // FIXME: double-check that we actually got it correct...
49
50 #ifdef _DEBUG
51 #ifdef WIN32
52 #undef THIS_FILE
53 static char THIS_FILE[] = __FILE__;
54 #endif
55 #endif
56
57 const char BSplineTensorSurface::ff_const_1[] = "BEGINBSPLINETENSORSURFACE";
58 const char BSplineTensorSurface::ff_const_2[] = "DIMENSIONU";
59 const char BSplineTensorSurface::ff_const_3[] = "DIMENSIONV";
60 const char BSplineTensorSurface::ff_const_4[] = "NUMBEROFCONTROLPOINTS";
61 const char BSplineTensorSurface::ff_const_5[] = "BEGINRATIONALBSPLINETENSORSURFACE";
62
63
64 int BSplineTensorSurface::CheckKnotPoints(const DCTPdvector& knots, int dim)
65 {
66 //now check knotvector whether it has
67 //(a,a,....a, ......., b,b,....b)
68 // |-dim+1-| |-dim+1-|
69 //structure
70
71 DCTPdvector::size_type max_index = knots.size() - 1;
72 double test_begin = knots[0], test_end = knots[max_index];
73
74 for(DCTPdvector::size_type i = 1; i < static_cast<unsigned int>(dim) + 1; ++i)
75 if(knots[i] != test_begin || knots[max_index - i] != test_end)
76 return -1; //FIXME: double comparison ?!
77
78 return 0;
79 }
80
81 int BSplineTensorSurface::deleteBezierKnot_U(double k)
82 {
83 DCTPdvector knots = basis_function_u.getKnotVector();
84
85 if(k >= knots[knots.size() - 1])
86 return -1; // knot is too high
87 if(k <= knots[0])
88 return -2; // knot is too low
89
90 DCTPdvector::size_type i = 0;
91 int mult = 0;
92 while(knots[i] <= k)
93 {
94 if(knots[i] == k)
95 mult++;
96 i++;
97 }
98 if(mult < dimension_u + 1)
99 return -3; // knot doesn't have enough multiplicity
100 i--;
101
102 // delete from knotvector
103 knots.erase(knots.begin() + i);
104 // set new knot vector
105 basis_function_u.setKnotVector(knots);
106 // delete control points from the control point mesh corresponding to just deleted knot
107 control_points.erase(control_points.begin() + i - dimension_u);
108 return 0;
109
110 }
111
112 int BSplineTensorSurface::deleteBezierKnot_V(double k)
113 {
114 DCTPdvector knots = basis_function_v.getKnotVector();
115
116 if(k >= knots[knots.size() - 1])
117 return -1; // knot is too high
118 if(k <= knots[0])
119 return -2; // knot is too low
120
121 DCTPdvector::size_type i = 0;
122 int mult = 0;
123 while(knots[i] <= k)
124 {
125 if(knots[i] == k)
126 mult++;
127 i++;
128 }
129 if(mult < dimension_v + 1)
130 return -3; // knot doesn't have enough multiplicity
131 i--;
132
133 // delete from knotvector
134 knots.erase(knots.begin() + i);
135 // set new knot vector
136 basis_function_v.setKnotVector(knots);
137
138 // delete control points from the control point mesh corresponding to just deleted knot
139 for(DCTPVec4dvector::size_type j = 0; j < control_points.size(); j++)
140 {
141 control_points[j].erase(control_points[j].begin() + i - dimension_v);
142 }
143
144 return 0;
145 }
146
147
148 //setup functions
149 int BSplineTensorSurface::setKnotsAndDimension(const DCTPdvector& knots_u, int dim_u, const DCTPdvector& knots_v, int dim_v)
150 {
151 if(dim_u < 1 || dim_v < 1)
152 return -1; //invalid dimension
153
154 DCTPdvector::size_type max_index = knots_u.size() - 1;
155 if(max_index < 3)
156 return -2; //here's an implicit check fer structure, see below
157 if(CheckKnotPoints(knots_u, dim_u) )
158 return -3;
159
160 max_index = knots_v.size() - 1;
161 if(max_index < 3)
162 return -2; //here's an implicit check fer structure, see below
163 if(CheckKnotPoints(knots_v, dim_v) )
164 return -3;
165
166 dimension_u = dim_u; dimension_v = dim_v;
167 if(basis_function_u.setKnotVector(knots_u) )
168 return -4; //error in BSplineBasisFunction.setKno...
169 if(basis_function_v.setKnotVector(knots_v) )
170 return -4; //error in BSplineBasisFunction.setKno...
171
172 return 0;
173 }
174
175 void BSplineTensorSurface::setControlPointMatrix(const DCTPVec4dmatrix &cps)
176 {
177 // FIXME: check that this is really a matrix
178 control_points = cps;
179 }
180
181 void BSplineTensorSurface::getParameterInterval_U(double &minpar, double &maxpar)
182 {
183 basis_function_u.getParameterInterval(minpar, maxpar);
184 }
185
186 void BSplineTensorSurface::getParameterInterval_V(double &minpar, double &maxpar)
187 {
188 basis_function_v.getParameterInterval(minpar, maxpar);
189 }
190
191 //I/O facilities - FIXME: read( char *fname ), etc. missing
192 //! Read the surface parameters from a given input stream.
193 /*!
194 *
195 * \param infile the stream from which the surface should be read.
196 * \return zero on success, a negative value on failure.
197 *
198 */
199 int BSplineTensorSurface::read(std::istream &infile)
200 {
201 //FIXME: maybe we need more checks!!!
202 char txtbuffer[256];
203 bool israt = false;
204
205
206 infile.getline(txtbuffer, 255); //read line
207 if(strcmp(txtbuffer, ff_const_1) &&
208 strcmp(txtbuffer, ff_const_5))
209 {
210 return -1; //bad file format
211 }
212 if(!strcmp(txtbuffer, ff_const_5))
213 {
214 israt = true;
215 }
216
217
218 infile >> txtbuffer; //FIXME: error prone: too long string causes problem!!!
219 if(strcmp(txtbuffer, ff_const_2) )
220 return -1; //yeah, bad file format again
221 infile >> dimension_u >> std::ws;
222 if(dimension_u < 1)
223 return -2; //ah, bad dimension
224
225 infile >> txtbuffer; //FIXME: error prone: too long string causes problem!!!
226 if(strcmp(txtbuffer, ff_const_3) )
227 return -1; //yeah, bad file format again
228 infile >> dimension_v >> std::ws;
229 if(dimension_v < 1)
230 return -2; //ah, bad dimension
231
232
233 if(basis_function_u.read(infile) )
234 return -3; //error reading basis function
235 if(CheckKnotPoints(basis_function_u.getKnotVector(), dimension_u) )
236 return -4;
237 infile >> std::ws;
238 if(basis_function_v.read(infile) )
239 return -3; //error reading basis function
240 if(CheckKnotPoints(basis_function_v.getKnotVector(), dimension_v) )
241 return -4;
242 infile >> std::ws;
243
244 infile >> txtbuffer; //FIXME: error prone: too long string causes problem!!!
245 if(strcmp(txtbuffer, ff_const_4) )
246 return -1; //bad file format once again
247
248 DCTPVec4dmatrix::size_type num_of_cps_u;
249 DCTPVec4dvector::size_type num_of_cps_v;
250 infile >> num_of_cps_u >> num_of_cps_v >> std::ws;
251 if(num_of_cps_u < 1 || num_of_cps_v < 1)
252 return -5; //too few control points
253 control_points.resize(num_of_cps_u); //FIXME: whatif not enoght memory?
254
255 for(DCTPdvector::size_type u = 0; u < num_of_cps_u; ++u)
256 {
257 control_points[u].resize(num_of_cps_v); //FIXME: not enough mem, exception thrown by STL
258
259 for(DCTPdvector::size_type v = 0; v < num_of_cps_v; ++v)
260 {
261 Vec4d cp;
262 if(israt)
263 {
264 infile >> cp[0] >> cp[1] >> cp[2] >> cp[3] >> std::ws;
265 }
266 else
267 {
268 infile >> cp[0] >> cp[1] >> cp[2] >> std::ws;
269 cp[3] = 1.0;
270 }
271 control_points[u][v] = cp; //FIXME: ya see, we need ERROR CHECKS!!!
272 }
273 }
274
275 return 0;
276 }
277
278 //! Write the surface parameters out to a given output stream.
279 /*!
280 *
281 * \param outfile the stream into which the surface should be written.
282 * \return zero on success, a negative value on failure.
283 *
284 */
285 int BSplineTensorSurface::write(std::ostream &outfile)
286 {
287 bool israt = false;
288 DCTPdvector::size_type u, v;
289 DCTPVec4dmatrix::size_type num_of_cps_u = control_points.size();
290 DCTPVec4dvector::size_type num_of_cps_v = control_points[0].size(); //FIXME:what if it's not initialized yet
291
292 for(u = 0; u < num_of_cps_u; ++u)
293 {
294 for(v = 0; v < num_of_cps_v; ++v)
295 {
296 if(osgAbs(control_points[u][v][3] - 1.0) > DCTP_EPS)
297 {
298 israt = true;
299 break;
300 }
301 }
302 }
303
304 //FIXME: maybe we need more checks!!!
305 outfile.precision(DCTP_PRECISION);
306 if(!israt)
307 {
308 outfile << ff_const_1 << std::endl;
309 }
310 else
311 {
312 outfile << ff_const_5 << std::endl;
313 }
314 outfile << ff_const_2 << " " << dimension_u << std::endl;
315 outfile << ff_const_3 << " " << dimension_v << std::endl;
316 basis_function_u.write(outfile);
317 basis_function_v.write(outfile);
318 outfile << ff_const_4 << " " << num_of_cps_u << " " << num_of_cps_v << std::endl;
319
320 for(u = 0; u < num_of_cps_u; ++u)
321 {
322 for(v = 0; v < num_of_cps_v; ++v)
323 {
324 Vec4d cp = control_points[u][v];
325 if(!israt)
326 {
327 outfile << cp[0] << " " << cp[1] << " " << cp[2] << std::endl;
328 }
329 else
330 {
331 outfile << cp[0] << " " << cp[1] << " " << cp[2] << " " << cp[3] << std::endl;
332 }
333 }
334 }
335
336 return 0;
337 }
338
339
340 //! Compute the surface at a given (u,v) point in parameter space.
341 /*!
342 *
343 * \param uv the point (in parameter space) where the surface should be
344 * computed.
345 * \param error flag to indicate if an error happened. It is set to 0
346 * if everything went fine, otherwise to a negative value.
347 * \return the value of the surface in 3D space.
348 *
349 */
350 Vec4d BSplineTensorSurface::compute4D(Vec2d uv, int &error)
351 {
352 //FIXME: verification before goin' into computation!!
353 double *nu;
354 int span_u = basis_function_u.computeAllNonzero(uv[0], dimension_u, nu);
355 double *nv;
356 int span_v = basis_function_v.computeAllNonzero(uv[1], dimension_v, nv);
357
358 Vec4d result(0.0, 0.0, 0.0, 0.0);
359 if(span_u < 0 || span_v < 0)
360 {
361 std::cerr << "spanu: " << span_u << " spanv: " << span_v << std::endl;
362 error = (span_u < 0) ? span_u : span_v;
363 if(span_u >= 0)
364 delete[] nu;
365 if(span_v >= 0)
366 delete[] nv;
367 result[3] = 1.0;
368 return result; //error occured in BSplineBasisFunction.computeAll...
369 }
370
371 unsigned int index_v = span_v - dimension_v;
372 Vec4d temp;
373 unsigned int index_u;
374
375 for(int l = 0; l <= dimension_v; ++l)
376 {
377 temp[0] = temp[1] = temp[2] = temp[3] = 0.0;
378 index_u = span_u - dimension_u;
379
380 for(DCTPVec4dmatrix::size_type k = 0; k <= static_cast<unsigned int>(dimension_u); ++k)
381 {
382 temp += control_points[index_u][index_v] * nu[k];
383 ++index_u;
384 }
385
386 ++index_v;
387 result += temp * nv[l];
388 }
389
390 delete[] nu;
391 delete[] nv;
392
393 return (result);
394 }
395
396 Vec3d BSplineTensorSurface::compute(Vec2d uv, int &error)
397 {
398 Vec4d rat_res;
399 Vec3d res;
400
401 rat_res = compute4D(uv, error);
402 res[0] = rat_res[0] / rat_res[3];
403 res[1] = rat_res[1] / rat_res[3];
404 res[2] = rat_res[2] / rat_res[3];
405 return res;
406 }
407
408 void BSplineTensorSurface::RatSurfaceDerivs(DCTPVec4dmatrix &rclHomDerivs, DCTPVec3dmatrix &rclEuclidDerivs)
409 {
410 unsigned int ui_k, ui_l, ui_j, ui_i;
411 Vec3d cl_v, cl_v2;
412
413 // FIXME: optimize for special case clHomDervis.size()==2!
414 // FIXME: binomial coefficients are all 1 in this case, so alredy removed
415
416 for(ui_k = 0; ui_k < rclHomDerivs.size(); ++ui_k)
417 {
418 for(ui_l = 0; ui_l < rclHomDerivs.size() - ui_k; ++ui_l)
419 {
420 cl_v[0] = rclHomDerivs[ui_k][ui_l][0];
421 cl_v[1] = rclHomDerivs[ui_k][ui_l][1];
422 cl_v[2] = rclHomDerivs[ui_k][ui_l][2];
423
424 for(ui_j = 1; ui_j <= ui_l; ++ui_j)
425 {
426 //cl_v -= Binomial(ui_l, ui_j) * rclHomDerivs[0][ui_j][3] * rclEuclidDerivs[ui_k][ui_l-ui_j];
427 cl_v -= rclHomDerivs[0][ui_j][3] * rclEuclidDerivs[ui_k][ui_l - ui_j];
428 }
429
430 for(ui_i = 1; ui_i <= ui_k; ++ui_i)
431 {
432 //cl_v -= Binomial(ui_k, ui_i) * rclHomDerivs[ui_i][0][3] * rclEuclidDerivs[ui_k-ui_i][ui_l];
433 cl_v -= rclHomDerivs[ui_i][0][3] * rclEuclidDerivs[ui_k - ui_i][ui_l];
434 cl_v2 = Vec3d(0.0, 0.0, 0.0);
435
436 for(ui_j = 1; ui_j <= ui_l; ++ui_j)
437 {
438 //cl_v2 += Binomial(ui_l, ui_j) * rclHomDerivs[ui_i][ui_j][3] * rclEuclidDerivs[ui_k-ui_i][ui_l-ui_j];
439 cl_v2 += rclHomDerivs[ui_i][ui_j][3] * rclEuclidDerivs[ui_k - ui_i][ui_l - ui_j];
440 }
441
442 // cl_v -= Binomial(ui_k, ui_i) * cl_v2;
443 cl_v -= cl_v2;
444 }
445
446 rclEuclidDerivs[ui_k][ui_l] = cl_v * (1.0 / rclHomDerivs[0][0][3]);
447 }
448 }
449 }
450
451
452 //! Compute the surface normal at a given (u,v) point in parameter space.
453 /*!
454 *
455 * \param uv the point (in parameter space) where the surface should be
456 * computed.
457 * \param error flag to indicate if an error happened. It is set to 0
458 * if everything went fine, otherwise to a negative value.
459 * \return the value of the surface in 3D space.
460 *
461 */
462 Vec3f BSplineTensorSurface::computeNormal(Vec2d clUV, int &riError, Pnt3f &rclPos)
463 {
464 Vec3d cl_norm;
465 Vec3d cl_pos;
466 Vec3f cl_normal;
467
468 cl_norm = computeNormal(clUV, riError, cl_pos);
469 cl_normal[0] = static_cast<Real32>(cl_norm[0]);
470 cl_normal[1] = static_cast<Real32>(cl_norm[1]);
471 cl_normal[2] = static_cast<Real32>(cl_norm[2]);
472 rclPos[0] = static_cast<Real32>(cl_pos[0]);
473 rclPos[1] = static_cast<Real32>(cl_pos[1]);
474 rclPos[2] = static_cast<Real32>(cl_pos[2]);
475
476 return cl_normal;
477 }
478
479
480 Vec3d BSplineTensorSurface::computeNormal(Vec2d clUV, int &riError, Vec3d &rclPos)
481 {
482 int i_uspan;
483 int i_vspan;
484 double **ppd_nu;
485 double **ppd_nv;
486 Vec3d cl_temp;
487 Vec3d cl_normal;
488 int i_k;
489 int i_s;
490 int i_l;
491 DCTPVec4dmatrix vvcl_skl;
492 DCTPVec3dmatrix vvcl_skl_eucl;
493 Vec4d *pcl_temp;
494 int i_r;
495 double d_len;
496 int i_u;
497 int i_v;
498
499 vvcl_skl.resize(2);
500 vvcl_skl[0].resize(2);
501 vvcl_skl[1].resize(2);
502 vvcl_skl_eucl.resize(2);
503 vvcl_skl_eucl[0].resize(2);
504 vvcl_skl_eucl[1].resize(2);
505
506 i_uspan = basis_function_u.computeDersBasisFuns(clUV[0], dimension_u, ppd_nu, 1);
507 i_vspan = basis_function_v.computeDersBasisFuns(clUV[1], dimension_v, ppd_nv, 1);
508 if( (i_uspan < 0) || (i_vspan < 0) )
509 {
510 riError = -1;
511 return cl_normal;
512 }
513
514 // vvcl_skl.resize( 2 );
515 // vvcl_skl[ 0 ].resize( 2 );
516 // vvcl_skl[ 1 ].resize( 2 );
517 pcl_temp = new Vec4d[dimension_v + 1];
518 // vcl_temp.resize( dimension_v + 1 );
519
520 for(i_k = 0; i_k <= 1; ++i_k)
521 {
522 i_v = i_vspan - dimension_v;
523
524 for(i_s = 0; i_s <= dimension_v; ++i_s)
525 {
526 pcl_temp[i_s] = Vec4d(0.0, 0.0, 0.0, 0.0);
527 i_u = i_uspan - dimension_u;
528
529 for(i_r = 0; i_r <= dimension_u; ++i_r)
530 {
531 pcl_temp[i_s] += control_points[i_u][i_v] * ppd_nu[i_k][i_r];
532 ++i_u;
533 }
534
535 ++i_v;
536 }
537
538 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
539 {
540 vvcl_skl[i_k][i_l] = Vec4d(0.0, 0.0, 0.0, 0.0);
541
542 for(i_s = 0; i_s <= dimension_v; ++i_s)
543 {
544 vvcl_skl[i_k][i_l] += pcl_temp[i_s] * ppd_nv[i_l][i_s];
545 }
546 }
547 }
548
549 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
550 correctDers(clUV, vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
551 cl_temp = vvcl_skl_eucl[1][0].cross(vvcl_skl_eucl[0][1]);
552 d_len = cl_temp.squareLength();
553 if(d_len < DCTP_EPS * DCTP_EPS)
554 {
555 // std::cerr<<"d_len too small: " << d_len << std::endl;
556 cl_normal[0] = cl_normal[1] = cl_normal[2] = 0.0;
557 riError = -2;
558 }
559 else
560 {
561 cl_temp *= 1.0 / sqrt(d_len);
562 riError = 0;
563 cl_normal[0] = cl_temp[0];
564 cl_normal[1] = cl_temp[1];
565 cl_normal[2] = cl_temp[2];
566 }
567
568 rclPos[0] = vvcl_skl_eucl[0][0][0];
569 rclPos[1] = vvcl_skl_eucl[0][0][1];
570 rclPos[2] = vvcl_skl_eucl[0][0][2];
571
572 // clean up allocated memory
573 delete[] ppd_nu[0];
574 delete[] ppd_nu[1];
575 delete[] ppd_nv[0];
576 delete[] ppd_nv[1];
577 delete[] ppd_nu;
578 delete[] ppd_nv;
579 delete[] pcl_temp;
580
581 return cl_normal;
582 }
583
584
585 //! Compute the surface normal and texture coordinate at a given (u,v) point in parameter space.
586 /*!
587 *
588 * \param uv the point (in parameter space) where the surface should be
589 * computed.
590 * \param error flag to indicate if an error happened. It is set to 0
591 * if everything went fine, otherwise to a negative value.
592 * \return the value of the surface in 3D space.
593 *
594 */
595 Vec3d BSplineTensorSurface::computeNormalTex(Vec2d & rclUV,
596 int & riError,
597 Vec3d & rclPos,
598 Vec2d & rclTexCoord,
599 const std::vector<std::vector<Pnt2f> > *cpvvclTexCP)
600 {
601 int i_uspan;
602 int i_vspan;
603 double **ppd_nu;
604 double **ppd_nv;
605 Vec3d cl_normal;
606 int i_k;
607 int i_s;
608 int i_l;
609 DCTPVec4dmatrix vvcl_skl;
610 DCTPVec3dmatrix vvcl_skl_eucl;
611 Vec4d *pcl_temp;
612 int i_r;
613 double d_len;
614 int i_u;
615 int i_v;
616 double *pd_nvl;
617 double *pd_nuk;
618 Vec2d cl_temp;
619
620 vvcl_skl.resize(2);
621 vvcl_skl[0].resize(2);
622 vvcl_skl[1].resize(2);
623 vvcl_skl_eucl.resize(2);
624 vvcl_skl_eucl[0].resize(2);
625 vvcl_skl_eucl[1].resize(2);
626
627 i_uspan = basis_function_u.computeDersBasisFuns(rclUV[0], dimension_u, ppd_nu, 1);
628 i_vspan = basis_function_v.computeDersBasisFuns(rclUV[1], dimension_v, ppd_nv, 1);
629 if( (i_uspan < 0) || (i_vspan < 0) )
630 {
631 riError = -1;
632 return cl_normal;
633 }
634
635 // vvcl_skl.resize( 2 );
636 // vvcl_skl[ 0 ].resize( 2 );
637 // vvcl_skl[ 1 ].resize( 2 );
638 pcl_temp = new Vec4d[dimension_v + 1];
639 // vcl_temp.resize( dimension_v + 1 );
640
641 for(i_k = 0; i_k <= 1; ++i_k)
642 {
643 i_v = i_vspan - dimension_v;
644
645 for(i_s = 0; i_s <= dimension_v; ++i_s)
646 {
647 pcl_temp[i_s] = Vec4d(0.0, 0.0, 0.0, 0.0);
648 i_u = i_uspan - dimension_u;
649
650 for(i_r = 0; i_r <= dimension_u; ++i_r)
651 {
652 pcl_temp[i_s] += control_points[i_u][i_v] * ppd_nu[i_k][i_r];
653 ++i_u;
654 }
655
656 ++i_v;
657 }
658
659 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
660 {
661 vvcl_skl[i_k][i_l] = Vec4d(0.0, 0.0, 0.0, 0.0);
662
663 for(i_s = 0; i_s <= dimension_v; ++i_s)
664 {
665 vvcl_skl[i_k][i_l] += pcl_temp[i_s] * ppd_nv[i_l][i_s];
666 }
667 }
668 }
669
670 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
671 rclPos = vvcl_skl_eucl[0][0];
672
673 correctDers(rclUV, vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
674 cl_normal = vvcl_skl_eucl[1][0].cross(vvcl_skl_eucl[0][1]);
675 d_len = cl_normal.squareLength();
676
677 if(d_len < DCTP_EPS * DCTP_EPS)
678 {
679 cl_normal[0] = cl_normal[1] = cl_normal[2] = 0.0;
680 riError = -2;
681 }
682 else
683 {
684 cl_normal *= 1.0 / sqrt(d_len);
685 riError = 0;
686 }
687
688 // calculate texture coord
689 rclTexCoord[0] = rclTexCoord[1] = 0.0;
690 i_v = i_vspan - dimension_v;
691 pd_nvl = ppd_nv[0];
692
693 for(i_l = 0; i_l <= dimension_v; ++i_l)
694 {
695 cl_temp[0] = cl_temp[1] = 0.0;
696 i_u = i_uspan - dimension_u;
697 pd_nuk = ppd_nu[0];
698
699 for(i_k = 0; i_k <= dimension_u; ++i_k)
700 {
701 cl_temp[0] += (*cpvvclTexCP)[i_u][i_v].x() * *pd_nuk;
702 cl_temp[1] += (*cpvvclTexCP)[i_u][i_v].y() * *pd_nuk;
703 ++i_u;
704 ++pd_nuk;
705 }
706
707 ++i_v;
708 rclTexCoord += cl_temp * *pd_nvl;
709 ++pd_nvl;
710 }
711
712 // clean up allocated memory
713 delete[] ppd_nu[0];
714 delete[] ppd_nu[1];
715 delete[] ppd_nv[0];
716 delete[] ppd_nv[1];
717 delete[] ppd_nu;
718 delete[] ppd_nv;
719 delete[] pcl_temp;
720
721 return cl_normal;
722 }
723
724
725 //! Insert a u knot into the surface. Returns < 0 on error, otherwise 0.
726 /*!
727 *
728 * \param k the knot to be inserted.
729 * \return zero on success, and a negative integer if some error occured.
730 *
731 */
732 int BSplineTensorSurface::insertKnot_U(double k)
733 {
734 DCTPdvector knots = basis_function_u.getKnotVector();
735 int span = basis_function_u.insertKnot(k);
736 if(span < 0)
737 return span; // some error happened
738
739 DCTPVec4dmatrix newcps(control_points.size() + 1);
740 DCTPVec4dvector::size_type i;
741
742 for(i = 0; i < newcps.size(); i++)
743 {
744 newcps[i].resize(control_points[0].size() ); // this must be 0 (or at least < newcps.size() - 1
745 }
746
747 for(i = 0; int(i) <= span - dimension_u; i++)
748 newcps[i] = control_points[i];
749
750 for(i = span - dimension_u + 1; i <= static_cast<unsigned int>(span); i++)
751 {
752 double alpha;
753 if(knots[i + dimension_u] != knots[i])
754 alpha = (k - knots[i]) / (knots[i + dimension_u] - knots[i]);
755 else
756 alpha = 0;
757
758 for(DCTPVec4dvector::size_type j = 0; j < newcps[i].size(); j++)
759 {
760 newcps[i][j] = control_points[i][j] * alpha +
761 control_points[i - 1][j] * (1 - alpha);
762 } // j
763
764 } // i
765
766 for(i = span + 1; i < newcps.size(); i++)
767 {
768 newcps[i] = control_points[i - 1];
769 }
770
771 control_points = newcps;
772 return span;
773 }
774
775
776 //! Insert v knot into the surface. Returns < 0 on error, otherwise 0.
777 /*!
778 *
779 * \param k the knot to be inserted.
780 * \return zero on success, and a negative integer if some error occured.
781 *
782 */
783 int BSplineTensorSurface::insertKnot_V(double k)
784 {
785 DCTPdvector knots = basis_function_v.getKnotVector();
786 int span = basis_function_v.insertKnot(k);
787 if(span < 0)
788 return span; // some error happened
789 DCTPVec4dmatrix newcps(control_points.size() );
790 DCTPVec4dvector::size_type i, j;
791
792 for(i = 0; i < newcps.size(); i++)
793 newcps[i].resize(control_points[i].size() + 1);
794
795 for(i = 0; i < control_points.size(); i++)
796 for(j = 0; int(j) <= span - dimension_v; j++)
797 newcps[i][j] = control_points[i][j];
798
799 // precalculate alpha values
800 DCTPdvector alphavec(dimension_v);
801
802 for(j = span - dimension_v + 1; int(j) <= span; j++)
803 {
804 if(knots[j + dimension_v] != knots[j])
805 alphavec[j - span + dimension_v - 1] = (k - knots[j]) / (knots[j + dimension_v] - knots[j]);
806 else
807 alphavec[j - span + dimension_v - 1] = 0;
808 }
809
810 for(i = 0; i < control_points.size(); i++)
811 for(j = span - dimension_v + 1; int(j) <= span; j++)
812 {
813 double alpha = alphavec[j - span + dimension_v - 1];
814 newcps[i][j] = control_points[i][j] * alpha + control_points[i][j - 1] * (1 - alpha);
815 }
816
817 // j
818
819 for(i = 0; i < control_points.size(); i++)
820 for(j = span + 1; j < newcps[i].size(); j++)
821 newcps[i][j] = control_points[i][j - 1];
822
823 control_points = newcps;
824 return span;
825 }
826
827 //! Convert a surface into Bezier patches.
828 /*!
829 * This function converts a surface into Bezier patches with
830 * appropriate knot insertions. (And if necessary knot deletions. )
831 * Leaves the original surface unchanged.
832 *
833 * \param beziers the matrix of BezierTensorSurface's to store the patches
834 * \param upars the u parameter vector to store where each patch begins in the
835 * original parameter space of the surface
836 * \param vpars the same for v parameter.
837 * \return zero on success, and a negative integer if some error occured.
838 *
839 * For example, for patch (i, j) it covers the parameter space of the
840 * original surface between upars[ i ]..upars[ i + 1 ] and
841 * vpars[ j ]..vpars[ j + 1 ].
842 */
843 int BSplineTensorSurface::makeBezier(beziersurfacematrix &beziers, DCTPdvector &upars, DCTPdvector &vpars)
844 {
845 double firstknotu, lastknotu;
846 double firstknotv, lastknotv;
847 basis_function_u.getParameterInterval(firstknotu, lastknotu);
848 basis_function_v.getParameterInterval(firstknotv, lastknotv);
849 DCTPdvector uknots = basis_function_u.getKnotVector();
850 DCTPdvector vknots = basis_function_v.getKnotVector();
851 DCTPdvector origuknots = uknots; // backup original curve characteristics
852 DCTPdvector origvknots = vknots; // same here
853 DCTPVec4dmatrix origcps = control_points; // same here
854 double prevknot = firstknotu;
855 int mult = 0;
856 int err;
857 int i;
858
859 // first convert along u knots
860 for(i = 1; i < int(uknots.size()); i++)
861 {
862 double actk = uknots[i];
863 if(actk == prevknot)
864 mult++;
865 else
866 {
867 for(int j = mult + 1; j < dimension_u; j++) //each interior knot must have the multiplicity of dimension -1
868 { //std::cerr << "about to insert uknot: " << prevknot << std::endl;
869 err = insertKnot_U(prevknot);
870 //std::cerr << "inserted uknot..." << std::endl;
871 if(err < 0)
872 return err; // some error happened during insertKnot
873 err = 0;
874 }
875
876 if(prevknot != firstknotu && prevknot != lastknotu)
877 {
878 for(int j = mult; j > dimension_u - 1; j--)
879 {
880 //std::cerr << "about to delete knot: " << prevknot << std::endl;
881 err = deleteBezierKnot_U(prevknot);
882 //std::cerr << "deleted knot..." << std::endl;
883 if(err)
884 return err; // some error happened during deleteBezierKnot
885 }
886 }
887 mult = 0;
888 }
889 prevknot = actk;
890 }
891
892 // now convert along v knots
893 mult = 0;
894 prevknot = firstknotv;
895
896 for(i = 1; i < int(vknots.size()); i++)
897 {
898 double actk = vknots[i];
899 if(actk == prevknot)
900 mult++;
901 else
902 {
903 for(int j = mult + 1; j < dimension_v; j++) //each interior knot must have the multiplicity of dimension -1
904 { // std::cerr << "about to insert vknot: " << prevknot << std::endl;
905 err = insertKnot_V(prevknot);
906 // std::cerr << "inserted vknot..." << std::endl;
907 if(err < 0)
908 return err; // some error happened during insertKnot
909 err = 0;
910 }
911
912 if(prevknot != firstknotv && prevknot != lastknotv)
913 {
914 for(int j = mult; j > dimension_v - 1; j--)
915 {
916 // std::cerr << "about to delete knot: " << prevknot << std::endl;
917 err = deleteBezierKnot_V(prevknot);
918 // std::cerr << "deleted knot..." << std::endl;
919 if(err)
920 return err; // some error happened during deleteBezierKnot
921 }
922 }
923 mult = 0;
924 }
925 prevknot = actk;
926 }
927
928 // now do the actual conversation into nxm bezier segments
929 uknots = basis_function_u.getKnotVector(); // FIXME: it could be more efficient.
930 int unum_of_beziers = (UInt32(uknots.size()) - 2) / dimension_u - 1;
931 if( (unum_of_beziers * dimension_u + 2) != int(uknots.size()) - dimension_u)
932 {
933 basis_function_u.setKnotVector(origuknots);
934 basis_function_v.setKnotVector(origuknots);
935 control_points = origcps;
936 return -5;
937 }
938 vknots = basis_function_v.getKnotVector(); // FIXME: it could be more efficient.
939 int vnum_of_beziers = (UInt32(vknots.size()) - 2) / dimension_v - 1;
940 if( (vnum_of_beziers * dimension_v + 2) != int(vknots.size()) - dimension_v)
941 {
942 basis_function_u.setKnotVector(origuknots);
943 basis_function_v.setKnotVector(origuknots);
944 control_points = origcps;
945 return -5;
946 }
947
948 beziers.resize(unum_of_beziers);
949 upars.resize(unum_of_beziers + 1);
950 vpars.resize(vnum_of_beziers + 1);
951
952 for(i = 0; i < unum_of_beziers + 1; i++)
953 upars[i] = uknots[1 + dimension_u * i];
954
955 for(i = 0; i < vnum_of_beziers + 1; i++)
956 vpars[i] = vknots[1 + dimension_v * i];
957
958 for(i = 0; i < unum_of_beziers; i++)
959 beziers[i].resize(vnum_of_beziers);
960
961 DCTPVec4dmatrix beziercps(dimension_u + 1);
962
963 for(i = 0; i < dimension_u + 1; i++)
964 beziercps[i].resize(dimension_v + 1);
965
966 for(i = 0; i < unum_of_beziers; i++)
967 {
968 for(int j = 0; j < vnum_of_beziers; j++)
969 {
970 for(int u = 0; u < dimension_u + 1; u++)
971 {
972 for(int v = 0; v < dimension_v + 1; v++)
973 {
974 beziercps[u][v] = control_points[i * dimension_u + u][j * dimension_v + v];
975 } // v
976
977 } // u
978
979 beziers[i][j].setControlPointMatrix(beziercps);
980 } // j
981
982 } // i
983
984 // restore original curve
985 basis_function_u.setKnotVector(origuknots);
986 basis_function_v.setKnotVector(origvknots);
987 control_points = origcps;
988 return 0;
989
990 }
991
992
993 int BSplineTensorSurface::getSplitParams(DCTPdvector &upars, DCTPdvector &vpars)
994 {
995 DCTPdvector &uknots = basis_function_u.getKnotVector();
996 DCTPdvector &vknots = basis_function_v.getKnotVector();
997 int unum_of_beziers = (UInt32(uknots.size()) - 2) / dimension_u - 1;
998 int vnum_of_beziers = (UInt32(vknots.size()) - 2) / dimension_v - 1;
999 int i;
1000
1001 upars.resize(unum_of_beziers + 1);
1002 vpars.resize(vnum_of_beziers + 1);
1003
1004 for(i = 0; i < unum_of_beziers + 1; i++)
1005 upars[i] = uknots[1 + dimension_u * i];
1006
1007 for(i = 0; i < vnum_of_beziers + 1; i++)
1008 vpars[i] = vknots[1 + dimension_v * i];
1009
1010 return 0;
1011 }
1012
1013
1014
1015 // subdivide surface at midpoint into 4 bezier surfaces
1016 int BSplineTensorSurface::midPointSubDivision(std::vector<std::vector<BSplineTensorSurface> > &rvvcl_newbspline)
1017 {
1018 double d_min;
1019 double d_max;
1020 int ui_cnt;
1021 int ui_idx;
1022 int i_span = 0;
1023 std::vector<std::vector<Vec4d> > vvcl_cp;
1024 unsigned int ui_cpy;
1025 unsigned int ui_cpy_cnt;
1026
1027 rvvcl_newbspline.resize(2);
1028 rvvcl_newbspline[0].resize(2);
1029 rvvcl_newbspline[1].resize(2);
1030
1031 // insert knots in u direction
1032 rvvcl_newbspline[0][0] = *this;
1033 rvvcl_newbspline[0][0].getParameterInterval_U(d_min, d_max);
1034 ui_cnt = rvvcl_newbspline[0][0].getDimension_U();
1035
1036 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1037 {
1038 i_span = rvvcl_newbspline[0][0].insertKnot_U( (d_min + d_max) * 0.5);
1039 }
1040
1041 i_span -= dimension_u - 2;
1042 rvvcl_newbspline[0][1].dimension_u = dimension_u;
1043 rvvcl_newbspline[1][0].dimension_u = dimension_u;
1044 rvvcl_newbspline[1][1].dimension_u = dimension_u;
1045 rvvcl_newbspline[0][1].dimension_v = dimension_v;
1046 rvvcl_newbspline[1][0].dimension_v = dimension_v;
1047 rvvcl_newbspline[1][1].dimension_v = dimension_v;
1048
1049 // split in u direction
1050 ui_cnt = UInt32(rvvcl_newbspline[0][0].getControlPointMatrix().size());
1051
1052 for(ui_idx = 0; ui_idx <= dimension_u; ++ui_idx)
1053 {
1054 rvvcl_newbspline[1][0].getKnotVector_U().push_back(
1055 rvvcl_newbspline[0][0].getKnotVector_U()[i_span]);
1056 }
1057
1058 for(ui_idx = i_span - 1; ui_idx < ui_cnt; ++ui_idx)
1059 {
1060 rvvcl_newbspline[1][0].getControlPointMatrix().push_back(
1061 rvvcl_newbspline[0][0].getControlPointMatrix()[ui_idx]);
1062 rvvcl_newbspline[1][0].getKnotVector_U().push_back(
1063 rvvcl_newbspline[0][0].getKnotVector_U()[ui_idx + dimension_u + 1]);
1064 }
1065
1066 rvvcl_newbspline[0][0].getControlPointMatrix().resize(i_span);
1067 rvvcl_newbspline[0][0].getKnotVector_U().resize(i_span + dimension_u);
1068 rvvcl_newbspline[0][0].getKnotVector_U().push_back(
1069 rvvcl_newbspline[0][0].getKnotVector_U()[i_span + dimension_u - 1]);
1070 rvvcl_newbspline[1][0].getKnotVector_V() = rvvcl_newbspline[0][0].getKnotVector_V();
1071
1072 ui_cnt = UInt32(rvvcl_newbspline[0][0].getKnotVector_U().size()) - 1;
1073 ui_idx = ui_cnt - dimension_u - 1;
1074 while(rvvcl_newbspline[0][0].getKnotVector_U()[ui_cnt] == rvvcl_newbspline[0][0].getKnotVector_U()[ui_idx])
1075 {
1076 rvvcl_newbspline[0][0].getKnotVector_U().pop_back();
1077 rvvcl_newbspline[0][0].getControlPointMatrix().pop_back();
1078 --ui_cnt;
1079 --ui_idx;
1080 }
1081
1082 // insert knots in v direction
1083 rvvcl_newbspline[0][0].getParameterInterval_V(d_min, d_max);
1084 ui_cnt = rvvcl_newbspline[0][0].getDimension_V();
1085
1086 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1087 {
1088 i_span = rvvcl_newbspline[0][0].insertKnot_V( (d_min + d_max) * 0.5);
1089 i_span = rvvcl_newbspline[1][0].insertKnot_V( (d_min + d_max) * 0.5);
1090 }
1091
1092 i_span -= dimension_v - 2;
1093
1094 // split in v direction
1095 ui_cnt = UInt32(rvvcl_newbspline[0][0].getControlPointMatrix()[0].size());
1096
1097 for(ui_idx = 0; ui_idx <= dimension_v; ++ui_idx)
1098 {
1099 rvvcl_newbspline[0][1].getKnotVector_V().push_back(
1100 rvvcl_newbspline[0][0].getKnotVector_V()[i_span]);
1101 rvvcl_newbspline[1][1].getKnotVector_V().push_back(
1102 rvvcl_newbspline[1][0].getKnotVector_V()[i_span]);
1103 }
1104
1105 for(ui_idx = i_span - 1; ui_idx < ui_cnt; ++ui_idx)
1106 {
1107 ui_cpy_cnt = UInt32(rvvcl_newbspline[0][0].getControlPointMatrix().size());
1108 rvvcl_newbspline[0][1].getControlPointMatrix().resize(ui_cpy_cnt);
1109
1110 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1111 {
1112 rvvcl_newbspline[0][1].getControlPointMatrix()[ui_cpy].push_back(
1113 rvvcl_newbspline[0][0].getControlPointMatrix()[ui_cpy][ui_idx]);
1114 }
1115
1116 ui_cpy_cnt = UInt32(rvvcl_newbspline[1][0].getControlPointMatrix().size());
1117 rvvcl_newbspline[1][1].getControlPointMatrix().resize(ui_cpy_cnt);
1118
1119 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1120 {
1121 rvvcl_newbspline[1][1].getControlPointMatrix()[ui_cpy].push_back(
1122 rvvcl_newbspline[1][0].getControlPointMatrix()[ui_cpy][ui_idx]);
1123 }
1124
1125 rvvcl_newbspline[0][1].getKnotVector_V().push_back(
1126 rvvcl_newbspline[0][0].getKnotVector_V()[ui_idx + dimension_v + 1]);
1127 rvvcl_newbspline[1][1].getKnotVector_V().push_back(
1128 rvvcl_newbspline[1][0].getKnotVector_V()[ui_idx + dimension_v + 1]);
1129 }
1130
1131 ui_cnt = UInt32(rvvcl_newbspline[0][0].getControlPointMatrix().size());
1132
1133 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1134 {
1135 rvvcl_newbspline[0][0].getControlPointMatrix()[ui_idx].resize(i_span);
1136 }
1137
1138 ui_cnt = UInt32(rvvcl_newbspline[1][0].getControlPointMatrix().size());
1139
1140 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1141 {
1142 rvvcl_newbspline[1][0].getControlPointMatrix()[ui_idx].resize(i_span);
1143 }
1144
1145 rvvcl_newbspline[0][0].getKnotVector_V().resize(i_span + dimension_v);
1146 rvvcl_newbspline[0][0].getKnotVector_V().push_back(
1147 rvvcl_newbspline[0][0].getKnotVector_V()[i_span + dimension_v - 1]);
1148 rvvcl_newbspline[1][0].getKnotVector_V().resize(i_span + dimension_v);
1149 rvvcl_newbspline[1][0].getKnotVector_V().push_back(
1150 rvvcl_newbspline[1][0].getKnotVector_V()[i_span + dimension_v - 1]);
1151 rvvcl_newbspline[0][1].getKnotVector_U() = rvvcl_newbspline[0][0].getKnotVector_U();
1152 rvvcl_newbspline[1][1].getKnotVector_U() = rvvcl_newbspline[1][0].getKnotVector_U();
1153
1154 ui_cnt = UInt32(rvvcl_newbspline[0][0].getKnotVector_V().size()) - 1;
1155 ui_idx = ui_cnt - dimension_v - 1;
1156 while(rvvcl_newbspline[0][0].getKnotVector_V()[ui_cnt] == rvvcl_newbspline[0][0].getKnotVector_V()[ui_idx])
1157 {
1158 rvvcl_newbspline[0][0].getKnotVector_V().pop_back();
1159 ui_cpy_cnt = UInt32(rvvcl_newbspline[0][0].getControlPointMatrix().size());
1160
1161 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1162 {
1163 rvvcl_newbspline[0][0].getControlPointMatrix()[ui_cpy].pop_back();
1164 }
1165
1166 rvvcl_newbspline[1][0].getKnotVector_V().pop_back();
1167 ui_cpy_cnt = UInt32(rvvcl_newbspline[1][0].getControlPointMatrix().size());
1168
1169 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1170 {
1171 rvvcl_newbspline[1][0].getControlPointMatrix()[ui_cpy].pop_back();
1172 }
1173
1174 --ui_cnt;
1175 --ui_idx;
1176 }
1177
1178 return 0;
1179 }
1180
1181
1182 // subdivide surface at midpoint into 2 bezier surfaces
1183 int BSplineTensorSurface::midPointSubDivisionU(std::vector<BSplineTensorSurface> &rvcl_newbspline)
1184 {
1185 double d_min;
1186 double d_max;
1187 int ui_cnt;
1188 int ui_idx;
1189 int i_span = 0;
1190 std::vector<std::vector<Vec4d> > vvcl_cp;
1191 // unsigned int ui_cpy;
1192 // unsigned int ui_cpy_cnt;
1193
1194 rvcl_newbspline.resize(2);
1195
1196 // insert knots in u direction
1197 rvcl_newbspline[0] = *this;
1198 getParameterInterval_U(d_min, d_max);
1199 ui_cnt = rvcl_newbspline[0].getDimension_U();
1200
1201 for(ui_idx = 0; ui_idx <= ui_cnt; ++ui_idx)
1202 {
1203 i_span = rvcl_newbspline[0].insertKnot_U( (d_min + d_max) * 0.5);
1204 }
1205
1206 i_span -= dimension_u - 2;
1207 rvcl_newbspline[1].dimension_u = dimension_u;
1208 rvcl_newbspline[1].dimension_v = dimension_v;
1209
1210 // split in u direction
1211 ui_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1212
1213 for(ui_idx = 0; ui_idx <= dimension_u; ++ui_idx)
1214 {
1215 rvcl_newbspline[1].getKnotVector_U().push_back(
1216 rvcl_newbspline[0].getKnotVector_U()[i_span]);
1217 }
1218
1219 for(ui_idx = i_span - 1; ui_idx < ui_cnt; ++ui_idx)
1220 {
1221 rvcl_newbspline[1].getControlPointMatrix().push_back(
1222 rvcl_newbspline[0].getControlPointMatrix()[ui_idx]);
1223 rvcl_newbspline[1].getKnotVector_U().push_back(
1224 rvcl_newbspline[0].getKnotVector_U()[ui_idx + dimension_u + 1]);
1225 }
1226
1227 rvcl_newbspline[0].getControlPointMatrix().resize(i_span);
1228 rvcl_newbspline[0].getKnotVector_U().resize(i_span + dimension_u);
1229 rvcl_newbspline[0].getKnotVector_U().push_back(
1230 rvcl_newbspline[0].getKnotVector_U()[i_span + dimension_u - 1]);
1231 rvcl_newbspline[1].getKnotVector_V() = rvcl_newbspline[0].getKnotVector_V();
1232
1233 ui_cnt = UInt32(rvcl_newbspline[0].getKnotVector_U().size()) - 1;
1234 ui_idx = ui_cnt - dimension_u - 1;
1235 while(rvcl_newbspline[0].getKnotVector_U()[ui_cnt] == rvcl_newbspline[0].getKnotVector_U()[ui_idx])
1236 {
1237 rvcl_newbspline[0].getKnotVector_U().pop_back();
1238 rvcl_newbspline[0].getControlPointMatrix().pop_back();
1239 --ui_cnt;
1240 --ui_idx;
1241 }
1242
1243 return 0;
1244 }
1245
1246
1247 // subdivide surface at midpoint into 2 bezier surfaces
1248 int BSplineTensorSurface::midPointSubDivisionV(std::vector<BSplineTensorSurface> &rvcl_newbspline)
1249 {
1250 double d_min;
1251 double d_max;
1252 int ui_cnt;
1253 int ui_idx;
1254 int i_span = 0;
1255 std::vector<std::vector<Vec4d> > vvcl_cp;
1256 unsigned int ui_cpy;
1257 unsigned int ui_cpy_cnt;
1258
1259 rvcl_newbspline.resize(2);
1260
1261 // insert knots in u direction
1262 rvcl_newbspline[0] = *this;
1263 getParameterInterval_V(d_min, d_max);
1264 ui_cnt = rvcl_newbspline[0].getDimension_V();
1265
1266 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1267 {
1268 i_span = rvcl_newbspline[0].insertKnot_V( (d_min + d_max) * 0.5);
1269 }
1270
1271 i_span -= dimension_v - 2;
1272 rvcl_newbspline[1].dimension_u = dimension_u;
1273 rvcl_newbspline[1].dimension_v = dimension_v;
1274
1275 // split in v direction
1276 ui_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix()[0].size());
1277
1278 for(ui_idx = 0; ui_idx <= dimension_v; ++ui_idx)
1279 {
1280 rvcl_newbspline[1].getKnotVector_V().push_back(
1281 rvcl_newbspline[0].getKnotVector_V()[i_span]);
1282 }
1283
1284 ui_cpy_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1285 rvcl_newbspline[1].getControlPointMatrix().resize(ui_cpy_cnt);
1286
1287 for(ui_idx = i_span - 1; ui_idx < ui_cnt; ++ui_idx)
1288 {
1289 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1290 {
1291 rvcl_newbspline[1].getControlPointMatrix()[ui_cpy].push_back(
1292 rvcl_newbspline[0].getControlPointMatrix()[ui_cpy][ui_idx]);
1293 }
1294
1295 rvcl_newbspline[1].getKnotVector_V().push_back(
1296 rvcl_newbspline[0].getKnotVector_V()[ui_idx + dimension_v + 1]);
1297 }
1298
1299 ui_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1300
1301 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1302 {
1303 rvcl_newbspline[0].getControlPointMatrix()[ui_idx].resize(i_span);
1304 }
1305
1306 rvcl_newbspline[0].getKnotVector_V().resize(i_span + dimension_v);
1307 rvcl_newbspline[0].getKnotVector_V().push_back(
1308 rvcl_newbspline[0].getKnotVector_V()[i_span + dimension_v - 1]);
1309 rvcl_newbspline[1].getKnotVector_U() = rvcl_newbspline[0].getKnotVector_U();
1310
1311 ui_cnt = UInt32(rvcl_newbspline[0].getKnotVector_V().size()) - 1;
1312 ui_idx = ui_cnt - dimension_v - 1;
1313 while(rvcl_newbspline[0].getKnotVector_V()[ui_cnt] == rvcl_newbspline[0].getKnotVector_V()[ui_idx])
1314 {
1315 rvcl_newbspline[0].getKnotVector_V().pop_back();
1316 ui_cpy_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1317
1318 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1319 {
1320 rvcl_newbspline[0].getControlPointMatrix()[ui_cpy].pop_back();
1321 }
1322
1323 --ui_cnt;
1324 --ui_idx;
1325 }
1326
1327 return 0;
1328 }
1329
1330
1331 // subdivide surface at param into 2 bezier surfaces
1332 int BSplineTensorSurface::subDivisionU(std::vector<BSplineTensorSurface> &rvcl_newbspline, double dSplit)
1333 {
1334 double d_min;
1335 double d_max;
1336 int ui_cnt;
1337 int ui_idx;
1338 int i_span = 0;
1339 std::vector<std::vector<Vec4d> > vvcl_cp;
1340 // unsigned int ui_cpy;
1341 // unsigned int ui_cpy_cnt;
1342
1343 rvcl_newbspline.clear();
1344 rvcl_newbspline.resize(2);
1345
1346 // insert knots in u direction
1347 rvcl_newbspline[0] = *this;
1348 getParameterInterval_U(d_min, d_max);
1349 ui_cnt = rvcl_newbspline[0].getDimension_U();
1350
1351 for(ui_idx = 0; ui_idx <= ui_cnt; ++ui_idx)
1352 {
1353 i_span = rvcl_newbspline[0].insertKnot_U(d_min + (d_max - d_min) * dSplit);
1354 }
1355
1356 i_span -= dimension_u - 2;
1357 rvcl_newbspline[1].dimension_u = dimension_u;
1358 rvcl_newbspline[1].dimension_v = dimension_v;
1359
1360 // split in u direction
1361 ui_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1362
1363 for(ui_idx = 0; ui_idx <= dimension_u; ++ui_idx)
1364 {
1365 rvcl_newbspline[1].getKnotVector_U().push_back(
1366 rvcl_newbspline[0].getKnotVector_U()[i_span]);
1367 }
1368
1369 for(ui_idx = i_span - 1; ui_idx < ui_cnt; ++ui_idx)
1370 {
1371 rvcl_newbspline[1].getControlPointMatrix().push_back(
1372 rvcl_newbspline[0].getControlPointMatrix()[ui_idx]);
1373 rvcl_newbspline[1].getKnotVector_U().push_back(
1374 rvcl_newbspline[0].getKnotVector_U()[ui_idx + dimension_u + 1]);
1375 }
1376
1377 rvcl_newbspline[0].getControlPointMatrix().resize(i_span);
1378 rvcl_newbspline[0].getKnotVector_U().resize(i_span + dimension_u);
1379 rvcl_newbspline[0].getKnotVector_U().push_back(
1380 rvcl_newbspline[0].getKnotVector_U()[i_span + dimension_u - 1]);
1381 rvcl_newbspline[1].getKnotVector_V() = rvcl_newbspline[0].getKnotVector_V();
1382
1383 ui_cnt = UInt32(rvcl_newbspline[0].getKnotVector_U().size()) - 1;
1384 ui_idx = ui_cnt - dimension_u - 1;
1385 while(rvcl_newbspline[0].getKnotVector_U()[ui_cnt] == rvcl_newbspline[0].getKnotVector_U()[ui_idx])
1386 {
1387 rvcl_newbspline[0].getKnotVector_U().pop_back();
1388 rvcl_newbspline[0].getControlPointMatrix().pop_back();
1389 --ui_cnt;
1390 --ui_idx;
1391 }
1392
1393 return 0;
1394 }
1395
1396
1397 // subdivide surface at param into 2 bezier surfaces
1398 int BSplineTensorSurface::subDivisionV(std::vector<BSplineTensorSurface> &rvcl_newbspline, double dSplit)
1399 {
1400 double d_min;
1401 double d_max;
1402 int ui_cnt;
1403 int ui_idx;
1404 int i_span = 0;
1405 std::vector<std::vector<Vec4d> > vvcl_cp;
1406 unsigned int ui_cpy;
1407 unsigned int ui_cpy_cnt;
1408
1409 rvcl_newbspline.clear();
1410 rvcl_newbspline.resize(2);
1411
1412 // insert knots in u direction
1413 rvcl_newbspline[0] = *this;
1414 getParameterInterval_V(d_min, d_max);
1415 ui_cnt = rvcl_newbspline[0].getDimension_V();
1416
1417 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1418 {
1419 i_span = rvcl_newbspline[0].insertKnot_V(d_min + (d_max - d_min) * dSplit);
1420 }
1421
1422 i_span -= dimension_v - 2;
1423 rvcl_newbspline[1].dimension_u = dimension_u;
1424 rvcl_newbspline[1].dimension_v = dimension_v;
1425
1426 // split in v direction
1427 ui_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix()[0].size());
1428
1429 for(ui_idx = 0; ui_idx <= dimension_v; ++ui_idx)
1430 {
1431 rvcl_newbspline[1].getKnotVector_V().push_back(
1432 rvcl_newbspline[0].getKnotVector_V()[i_span]);
1433 }
1434
1435 ui_cpy_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1436 rvcl_newbspline[1].getControlPointMatrix().resize(ui_cpy_cnt);
1437
1438 for(ui_idx = i_span - 1; ui_idx < ui_cnt; ++ui_idx)
1439 {
1440 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1441 {
1442 rvcl_newbspline[1].getControlPointMatrix()[ui_cpy].push_back(
1443 rvcl_newbspline[0].getControlPointMatrix()[ui_cpy][ui_idx]);
1444 }
1445
1446 rvcl_newbspline[1].getKnotVector_V().push_back(
1447 rvcl_newbspline[0].getKnotVector_V()[ui_idx + dimension_v + 1]);
1448 }
1449
1450 ui_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1451
1452 for(ui_idx = 0; ui_idx < ui_cnt; ++ui_idx)
1453 {
1454 rvcl_newbspline[0].getControlPointMatrix()[ui_idx].resize(i_span);
1455 }
1456
1457 rvcl_newbspline[0].getKnotVector_V().resize(i_span + dimension_v);
1458 rvcl_newbspline[0].getKnotVector_V().push_back(
1459 rvcl_newbspline[0].getKnotVector_V()[i_span + dimension_v - 1]);
1460 rvcl_newbspline[1].getKnotVector_U() = rvcl_newbspline[0].getKnotVector_U();
1461
1462 ui_cnt = UInt32(rvcl_newbspline[0].getKnotVector_V().size()) - 1;
1463 ui_idx = ui_cnt - dimension_v - 1;
1464 while(rvcl_newbspline[0].getKnotVector_V()[ui_cnt] == rvcl_newbspline[0].getKnotVector_V()[ui_idx])
1465 {
1466 rvcl_newbspline[0].getKnotVector_V().pop_back();
1467 ui_cpy_cnt = UInt32(rvcl_newbspline[0].getControlPointMatrix().size());
1468
1469 for(ui_cpy = 0; ui_cpy < ui_cpy_cnt; ++ui_cpy)
1470 {
1471 rvcl_newbspline[0].getControlPointMatrix()[ui_cpy].pop_back();
1472 }
1473
1474 --ui_cnt;
1475 --ui_idx;
1476 }
1477
1478 return 0;
1479 }
1480
1481
1482 //new vector-enabled computational functions from Michael
1483 void BSplineTensorSurface::compute(std::vector<Vec2d> &rvclUV,
1484 std::vector<Pnt3f> &rvclPos)
1485
1486 {
1487 const unsigned int cui_cnt = UInt32(rvclUV.size());
1488 unsigned int ui_idx;
1489 double *pd_nu;
1490 double *pd_nv;
1491 int i_span_u = -1;
1492 int i_span_v = -1;
1493 unsigned int ui_index_v;
1494 Vec4d *pcl_temp = new Vec4d[dimension_u + 1];
1495 unsigned int ui_index_u;
1496 int i_l;
1497 int i_k;
1498 double *pd_nul;
1499 double *pd_nvk;
1500 bool b_u_new;
1501 bool b_v_new;
1502 Vec4d *pcl_templ;
1503 int i_old_u;
1504
1505 rvclPos.resize(cui_cnt);
1506
1507 for(ui_idx = 0; ui_idx < cui_cnt; ++ui_idx)
1508 {
1509 if(ui_idx == 0)
1510 {
1511 i_span_u = basis_function_u.computeAllNonzero(rvclUV[ui_idx][0], dimension_u, pd_nu);
1512 i_span_v = basis_function_v.computeAllNonzero(rvclUV[ui_idx][1], dimension_v, pd_nv);
1513 b_u_new = b_v_new = true;
1514 }
1515 else
1516 {
1517 if(osgAbs(rvclUV[ui_idx][0] - rvclUV[ui_idx - 1][0]) > DCTP_EPS)
1518 {
1519 i_old_u = i_span_u;
1520 if(i_span_u >= 0)
1521 delete[] pd_nu;
1522 i_span_u = basis_function_u.computeAllNonzero(rvclUV[ui_idx][0], dimension_u, pd_nu);
1523 b_u_new = true;
1524 b_v_new = (i_old_u != i_span_u) ? true : false;
1525 }
1526 else
1527 {
1528 b_u_new = false;
1529 b_v_new = false;
1530 }
1531 if(osgAbs(rvclUV[ui_idx][1] - rvclUV[ui_idx - 1][1]) > DCTP_EPS)
1532 {
1533 if(i_span_v >= 0)
1534 delete[] pd_nv;
1535 i_span_v = basis_function_v.computeAllNonzero(rvclUV[ui_idx][1], dimension_v, pd_nv);
1536 b_v_new = true;
1537 }
1538 }
1539
1540 if( (i_span_u < 0) || (i_span_v < 0) )
1541 {
1542 rvclPos[ui_idx] = Pnt3f(0.0, 0.0, 0.0);
1543 }
1544 else if(b_v_new)
1545 {
1546 ui_index_u = i_span_u - dimension_u;
1547 pd_nul = pd_nu;
1548 pcl_templ = pcl_temp;
1549 Vec4d tempposition(0.0, 0.0, 0.0, 0.0);
1550
1551 for(i_l = 0; i_l <= dimension_u; ++i_l)
1552 {
1553 (*pcl_templ) = Vec4d(0.0, 0.0, 0.0, 0.0);
1554 ui_index_v = i_span_v - dimension_v;
1555 pd_nvk = pd_nv;
1556
1557 for(i_k = 0; i_k <= dimension_v; ++i_k)
1558 {
1559 *pcl_templ += control_points[ui_index_u][ui_index_v] * *pd_nvk;
1560 ++ui_index_v;
1561 ++pd_nvk;
1562 }
1563
1564 ++ui_index_u;
1565 tempposition += *pcl_templ * *pd_nul;
1566 ++pd_nul;
1567 ++pcl_templ;
1568 }
1569
1570 rvclPos[ui_idx][0] = tempposition[0] / tempposition[3];
1571 rvclPos[ui_idx][1] = tempposition[1] / tempposition[3];
1572 rvclPos[ui_idx][2] = tempposition[2] / tempposition[3];
1573 }
1574 else if(b_u_new)
1575 {
1576 ui_index_u = i_span_u - dimension_u;
1577 pd_nul = pd_nu;
1578 pcl_templ = pcl_temp;
1579 Vec4d tempposition(0.0, 0.0, 0.0, 0.0);
1580
1581 for(i_l = 0; i_l <= dimension_u; ++i_l)
1582 {
1583 ++ui_index_u;
1584 tempposition += *pcl_templ * *pd_nul;
1585 ++pd_nul;
1586 ++pcl_templ;
1587 }
1588
1589 rvclPos[ui_idx][0] = tempposition[0] / tempposition[3];
1590 rvclPos[ui_idx][1] = tempposition[1] / tempposition[3];
1591 rvclPos[ui_idx][2] = tempposition[2] / tempposition[3];
1592 #ifdef OSG_COUNT_COMP
1593 ++g_uiVSameComp;
1594 #endif
1595 }
1596 else
1597 {
1598 rvclPos[ui_idx] = rvclPos[ui_idx - 1];
1599 #ifdef OSG_COUNT_COMP
1600 ++g_uiSameComp;
1601 #endif
1602 }
1603 #ifdef OSG_COUNT_COMP
1604 ++g_uiTotalComp;
1605 #endif
1606 }
1607
1608 if(i_span_u >= 0)
1609 delete[] pd_nu;
1610 if(i_span_v >= 0)
1611 delete[] pd_nv;
1612 delete[] pcl_temp;
1613
1614 #ifdef OSG_COUNT_COMP
1615 cerr << g_uiTotalComp << " " << g_uiVSameComp << " " << g_uiSameComp << endl;
1616 #endif
1617 }
1618
1619
1620 void BSplineTensorSurface::computeNormal(std::vector<Vec2d> &rvclUV,
1621 std::vector<Pnt3f> &rvclPos,
1622 std::vector<Vec3f> &rvclNorm)
1623 {
1624 const unsigned int cui_cnt = UInt32(rvclUV.size());
1625 unsigned int ui_idx;
1626 int i_uspan = -1;
1627 int i_vspan = -1;
1628 double **ppd_nu;
1629 double **ppd_nv;
1630 int i_k;
1631 int i_s;
1632 int i_l;
1633 DCTPVec4dmatrix vvcl_skl;
1634 DCTPVec3dmatrix vvcl_skl_eucl;
1635 Vec4d *apcl_temp[2];
1636 int i_r;
1637 double d_len = 0.0;
1638 int i_u;
1639 int i_v;
1640 bool b_u_new;
1641 bool b_v_new;
1642 Vec4d *pcl_temps;
1643 double *pd_nuls;
1644 double *pd_nvkr;
1645 Vec4d *pcl_sklkl;
1646 int i_old_u;
1647
1648 // std::cerr << dimension_u << " " << dimension_v << std::endl;
1649
1650 vvcl_skl.resize(2);
1651 vvcl_skl[0].resize(2);
1652 vvcl_skl[1].resize(2);
1653 vvcl_skl_eucl.resize(2);
1654 vvcl_skl_eucl[0].resize(2);
1655 vvcl_skl_eucl[1].resize(2);
1656
1657 apcl_temp[0] = new Vec4d[dimension_u + 1];
1658 apcl_temp[1] = new Vec4d[dimension_u + 1];
1659
1660 rvclPos.resize(cui_cnt);
1661 rvclNorm.resize(cui_cnt);
1662
1663 for(ui_idx = 0; ui_idx < cui_cnt; ++ui_idx)
1664 {
1665 if(ui_idx == 0)
1666 {
1667 i_uspan = basis_function_u.computeDersBasisFuns(rvclUV[ui_idx][0], dimension_u, ppd_nu, 1);
1668 i_vspan = basis_function_v.computeDersBasisFuns(rvclUV[ui_idx][1], dimension_v, ppd_nv, 1);
1669 b_u_new = b_v_new = true;
1670 }
1671 else
1672 {
1673 if(osgAbs(rvclUV[ui_idx][0] - rvclUV[ui_idx - 1][0]) > DCTP_EPS)
1674 {
1675 i_old_u = i_uspan;
1676 if(i_uspan >= 0)
1677 {
1678 delete[] ppd_nu[0];
1679 delete[] ppd_nu[1];
1680 delete[] ppd_nu;
1681 }
1682 i_uspan = basis_function_u.computeDersBasisFuns(rvclUV[ui_idx][0], dimension_u, ppd_nu, 1);
1683 b_u_new = true;
1684 b_v_new = (i_old_u != i_uspan) ? true : false;
1685 }
1686 else
1687 {
1688 b_u_new = false;
1689 b_v_new = false;
1690 }
1691 if(osgAbs(rvclUV[ui_idx][1] - rvclUV[ui_idx - 1][1]) > DCTP_EPS)
1692 {
1693 if(i_vspan >= 0)
1694 {
1695 delete[] ppd_nv[0];
1696 delete[] ppd_nv[1];
1697 delete[] ppd_nv;
1698 }
1699 i_vspan = basis_function_v.computeDersBasisFuns(rvclUV[ui_idx][1], dimension_v, ppd_nv, 1);
1700 b_v_new = true;
1701 }
1702 }
1703 if( (i_uspan < 0) || (i_vspan < 0) )
1704 {
1705 rvclPos[ui_idx] = Pnt3f(0.0, 0.0, 0.0);
1706 rvclNorm[ui_idx] = Vec3f(0.0, 0.0, 0.0);
1707 }
1708 else if(b_v_new)
1709 {
1710 for(i_k = 0; i_k <= 1; ++i_k)
1711 {
1712 i_u = i_uspan - dimension_u;
1713 pcl_temps = apcl_temp[i_k];
1714
1715 for(i_s = 0; i_s <= dimension_u; ++i_s)
1716 {
1717 (*pcl_temps) = Vec4d(0.0, 0.0, 0.0, 0.0);
1718 i_v = i_vspan - dimension_v;
1719 pd_nvkr = ppd_nv[i_k];
1720
1721 for(i_r = 0; i_r <= dimension_v; ++i_r)
1722 {
1723 *pcl_temps += control_points[i_u][i_v] * *pd_nvkr;
1724 ++i_v;
1725 ++pd_nvkr;
1726 }
1727
1728 ++i_u;
1729 ++pcl_temps;
1730 }
1731
1732 pcl_sklkl = &(vvcl_skl[i_k][0]);
1733
1734 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
1735 {
1736 (*pcl_sklkl) = Vec4d(0.0, 0.0, 0.0, 0.0);
1737 pcl_temps = apcl_temp[i_k];
1738 pd_nuls = ppd_nu[i_l];
1739
1740 for(i_s = 0; i_s <= dimension_u; ++i_s)
1741 {
1742 *pcl_sklkl += *pcl_temps * *pd_nuls;
1743 ++pcl_temps;
1744 ++pd_nuls;
1745 }
1746
1747 ++pcl_sklkl;
1748 }
1749 }
1750
1751 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
1752 rvclPos[ui_idx][0] = vvcl_skl_eucl[0][0][0];
1753 rvclPos[ui_idx][1] = vvcl_skl_eucl[0][0][1];
1754 rvclPos[ui_idx][2] = vvcl_skl_eucl[0][0][2];
1755
1756 correctDers(rvclUV[ui_idx], vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
1757 // rvclNorm[ ui_idx ] = aacl_skl[ 0 ][ 1 ].cross( aacl_skl[ 1 ][ 0 ] ); // switched because of optimization!
1758 // d_len = rvclNorm[ ui_idx ].squareLength( );
1759 Vec3d tempnorm = vvcl_skl_eucl[0][1].cross(vvcl_skl_eucl[1][0]);
1760 rvclNorm[ui_idx][0] = tempnorm[0];
1761 rvclNorm[ui_idx][1] = tempnorm[1];
1762 rvclNorm[ui_idx][2] = tempnorm[2];
1763 d_len = tempnorm.squareLength();
1764
1765 if(d_len > DCTP_EPS * DCTP_EPS)
1766 {
1767 rvclNorm[ui_idx] *= 1.0 / sqrt(d_len);
1768 }
1769 else
1770 {
1771 rvclNorm[ui_idx][0] = rvclNorm[ui_idx][1] = rvclNorm[ui_idx][2] = 0.0;
1772 }
1773 }
1774 else if(b_u_new)
1775 {
1776 for(i_k = 0; i_k <= 1; ++i_k)
1777 {
1778 pcl_sklkl = &(vvcl_skl[i_k][0]);
1779
1780 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
1781 {
1782 (*pcl_sklkl) = Vec4d(0.0, 0.0, 0.0, 0.0);
1783 pcl_temps = apcl_temp[i_k];
1784 pd_nuls = ppd_nu[i_l];
1785
1786 for(i_s = 0; i_s <= dimension_u; ++i_s)
1787 {
1788 *pcl_sklkl += *pcl_temps * *pd_nuls;
1789 ++pcl_temps;
1790 ++pd_nuls;
1791 }
1792
1793 ++pcl_sklkl;
1794 }
1795 }
1796
1797 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
1798 rvclPos[ui_idx][0] = vvcl_skl_eucl[0][0][0];
1799 rvclPos[ui_idx][1] = vvcl_skl_eucl[0][0][1];
1800 rvclPos[ui_idx][2] = vvcl_skl_eucl[0][0][2];
1801
1802 correctDers(rvclUV[ui_idx], vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
1803 // rvclNorm[ ui_idx ] = aacl_skl[ 0 ][ 1 ].cross( aacl_skl[ 1 ][ 0 ] ); // switched because of optimization!
1804 // d_len = rvclNorm[ ui_idx ].squareLength( );
1805 Vec3d tempnorm = vvcl_skl_eucl[0][1].cross(vvcl_skl_eucl[1][0]);
1806 rvclNorm[ui_idx][0] = tempnorm[0];
1807 rvclNorm[ui_idx][1] = tempnorm[1];
1808 rvclNorm[ui_idx][2] = tempnorm[2];
1809 if(d_len > DCTP_EPS * DCTP_EPS)
1810 {
1811 rvclNorm[ui_idx] *= 1.0 / sqrt(d_len);
1812 }
1813 else
1814 {
1815 rvclNorm[ui_idx][0] = rvclNorm[ui_idx][1] = rvclNorm[ui_idx][2] = 0.0;
1816 }
1817 #ifdef OSG_COUNT_COMP
1818 ++g_uiVSameComp;
1819 #endif
1820 }
1821 else
1822 {
1823 rvclPos[ui_idx] = rvclPos[ui_idx - 1];
1824 rvclNorm[ui_idx] = rvclNorm[ui_idx - 1];
1825 #ifdef OSG_COUNT_COMP
1826 ++g_uiSameComp;
1827 #endif
1828 }
1829 #ifdef OSG_COUNT_COMP
1830 ++g_uiTotalComp;
1831 #endif
1832 }
1833
1834 // clean up allocated memory
1835 if(i_uspan >= 0)
1836 {
1837 delete[] ppd_nu[0];
1838 delete[] ppd_nu[1];
1839 delete[] ppd_nu;
1840 }
1841 if(i_vspan >= 0)
1842 {
1843 delete[] ppd_nv[0];
1844 delete[] ppd_nv[1];
1845 delete[] ppd_nv;
1846 }
1847 delete[] apcl_temp[0];
1848 delete[] apcl_temp[1];
1849
1850 #ifdef OSG_COUNT_COMP
1851 cerr << g_uiTotalComp << " " << g_uiVSameComp << " " << g_uiSameComp << endl;
1852 #endif
1853 }
1854
1855
1856 void BSplineTensorSurface::computeNormalforTrimming(std::vector<Vec2d> &rvclUV,
1857 std::vector<Vec3d> &rvclNorm,
1858 std::vector<Vec3d> *pvclPos)
1859 {
1860 const unsigned int cui_cnt = UInt32(rvclUV.size());
1861 unsigned int ui_idx;
1862 int i_uspan = -1;
1863 int i_vspan = -1;
1864 double **ppd_nu;
1865 double **ppd_nv;
1866 int i_k;
1867 int i_s;
1868 int i_l;
1869 DCTPVec4dmatrix vvcl_skl;
1870 DCTPVec3dmatrix vvcl_skl_eucl;
1871 Vec4d *apcl_temp[2];
1872 int i_r;
1873 double d_len = 0.0;
1874 int i_u;
1875 int i_v;
1876 bool b_u_new;
1877 bool b_v_new;
1878 Vec4d *pcl_temps;
1879 double *pd_nuls;
1880 double *pd_nvkr;
1881 Vec4d *pcl_sklkl;
1882 int i_old_u;
1883
1884 // std::cerr << dimension_u << " " << dimension_v << std::endl;
1885
1886 vvcl_skl.resize(2);
1887 vvcl_skl[0].resize(2);
1888 vvcl_skl[1].resize(2);
1889 vvcl_skl_eucl.resize(2);
1890 vvcl_skl_eucl[0].resize(2);
1891 vvcl_skl_eucl[1].resize(2);
1892
1893 apcl_temp[0] = new Vec4d[dimension_u + 1];
1894 apcl_temp[1] = new Vec4d[dimension_u + 1];
1895
1896 if(pvclPos != NULL)
1897 {
1898 pvclPos->resize(cui_cnt);
1899 }
1900 rvclNorm.resize(cui_cnt);
1901
1902 for(ui_idx = 0; ui_idx < cui_cnt; ++ui_idx)
1903 {
1904 if(ui_idx == 0)
1905 {
1906 i_uspan = basis_function_u.computeDersBasisFuns(rvclUV[ui_idx][0], dimension_u, ppd_nu, 1);
1907 i_vspan = basis_function_v.computeDersBasisFuns(rvclUV[ui_idx][1], dimension_v, ppd_nv, 1);
1908 b_u_new = b_v_new = true;
1909 }
1910 else
1911 {
1912 if(osgAbs(rvclUV[ui_idx][0] - rvclUV[ui_idx - 1][0]) > DCTP_EPS)
1913 {
1914 i_old_u = i_uspan;
1915 if(i_uspan >= 0)
1916 {
1917 delete[] ppd_nu[0];
1918 delete[] ppd_nu[1];
1919 delete[] ppd_nu;
1920 }
1921 i_uspan = basis_function_u.computeDersBasisFuns(rvclUV[ui_idx][0], dimension_u, ppd_nu, 1);
1922 b_u_new = true;
1923 b_v_new = (i_old_u != i_uspan) ? true : false;
1924 }
1925 else
1926 {
1927 b_u_new = false;
1928 b_v_new = false;
1929 }
1930 if(osgAbs(rvclUV[ui_idx][1] - rvclUV[ui_idx - 1][1]) > DCTP_EPS)
1931 {
1932 if(i_vspan >= 0)
1933 {
1934 delete[] ppd_nv[0];
1935 delete[] ppd_nv[1];
1936 delete[] ppd_nv;
1937 }
1938 i_vspan = basis_function_v.computeDersBasisFuns(rvclUV[ui_idx][1], dimension_v, ppd_nv, 1);
1939 b_v_new = true;
1940 }
1941 }
1942 if( (i_uspan < 0) || (i_vspan < 0) )
1943 {
1944 if(pvclPos != NULL)
1945 {
1946 (*pvclPos)[ui_idx] = Vec3d(0.0, 0.0, 0.0);
1947 }
1948 rvclNorm[ui_idx] = Vec3d(0.0, 0.0, 0.0);
1949 }
1950 else if(b_v_new)
1951 {
1952 for(i_k = 0; i_k <= 1; ++i_k)
1953 {
1954 i_u = i_uspan - dimension_u;
1955 pcl_temps = apcl_temp[i_k];
1956
1957 for(i_s = 0; i_s <= dimension_u; ++i_s)
1958 {
1959 (*pcl_temps) = Vec4d(0.0, 0.0, 0.0, 0.0);
1960 i_v = i_vspan - dimension_v;
1961 pd_nvkr = ppd_nv[i_k];
1962
1963 for(i_r = 0; i_r <= dimension_v; ++i_r)
1964 {
1965 *pcl_temps += control_points[i_u][i_v] * *pd_nvkr;
1966 ++i_v;
1967 ++pd_nvkr;
1968 }
1969
1970 ++i_u;
1971 ++pcl_temps;
1972 }
1973
1974 pcl_sklkl = &(vvcl_skl[i_k][0]);
1975
1976 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
1977 {
1978 (*pcl_sklkl) = Vec4d(0.0, 0.0, 0.0, 0.0);
1979 pcl_temps = apcl_temp[i_k];
1980 pd_nuls = ppd_nu[i_l];
1981
1982 for(i_s = 0; i_s <= dimension_u; ++i_s)
1983 {
1984 *pcl_sklkl += *pcl_temps * *pd_nuls;
1985 ++pcl_temps;
1986 ++pd_nuls;
1987 }
1988
1989 ++pcl_sklkl;
1990 }
1991 }
1992
1993 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
1994 if(pvclPos != NULL)
1995 {
1996 (*pvclPos)[ui_idx] = vvcl_skl_eucl[0][0];
1997 }
1998
1999 correctDers(rvclUV[ui_idx], vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
2000 // rvclNorm[ ui_idx ] = aacl_skl[ 0 ][ 1 ].cross( aacl_skl[ 1 ][ 0 ] ); // switched because of optimization!
2001 // d_len = rvclNorm[ ui_idx ].squareLength( );
2002 Vec3d tempnorm = vvcl_skl_eucl[0][1].cross(vvcl_skl_eucl[1][0]);
2003 rvclNorm[ui_idx][0] = tempnorm[0];
2004 rvclNorm[ui_idx][1] = tempnorm[1];
2005 rvclNorm[ui_idx][2] = tempnorm[2];
2006 d_len = tempnorm.squareLength();
2007
2008 if(d_len > DCTP_EPS * DCTP_EPS)
2009 {
2010 rvclNorm[ui_idx] *= 1.0 / sqrt(d_len);
2011 }
2012 else
2013 {
2014 rvclNorm[ui_idx][0] = rvclNorm[ui_idx][1] = rvclNorm[ui_idx][2] = 0.0;
2015 }
2016 }
2017 else if(b_u_new)
2018 {
2019 for(i_k = 0; i_k <= 1; ++i_k)
2020 {
2021 pcl_sklkl = &(vvcl_skl[i_k][0]);
2022
2023 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
2024 {
2025 (*pcl_sklkl) = Vec4d(0.0, 0.0, 0.0, 0.0);
2026 pcl_temps = apcl_temp[i_k];
2027 pd_nuls = ppd_nu[i_l];
2028
2029 for(i_s = 0; i_s <= dimension_u; ++i_s)
2030 {
2031 *pcl_sklkl += *pcl_temps * *pd_nuls;
2032 ++pcl_temps;
2033 ++pd_nuls;
2034 }
2035
2036 ++pcl_sklkl;
2037 }
2038 }
2039
2040 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
2041 if(pvclPos != NULL)
2042 {
2043 (*pvclPos)[ui_idx] = vvcl_skl_eucl[0][0];
2044 }
2045
2046 correctDers(rvclUV[ui_idx], vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
2047 // rvclNorm[ ui_idx ] = aacl_skl[ 0 ][ 1 ].cross( aacl_skl[ 1 ][ 0 ] ); // switched because of optimization!
2048 // d_len = rvclNorm[ ui_idx ].squareLength( );
2049 Vec3d tempnorm = vvcl_skl_eucl[0][1].cross(vvcl_skl_eucl[1][0]);
2050 rvclNorm[ui_idx][0] = tempnorm[0];
2051 rvclNorm[ui_idx][1] = tempnorm[1];
2052 rvclNorm[ui_idx][2] = tempnorm[2];
2053 if(d_len > DCTP_EPS * DCTP_EPS)
2054 {
2055 rvclNorm[ui_idx] *= 1.0 / sqrt(d_len);
2056 }
2057 else
2058 {
2059 rvclNorm[ui_idx][0] = rvclNorm[ui_idx][1] = rvclNorm[ui_idx][2] = 0.0;
2060 }
2061 #ifdef OSG_COUNT_COMP
2062 ++g_uiVSameComp;
2063 #endif
2064 }
2065 else
2066 {
2067 if(pvclPos != NULL)
2068 {
2069 (*pvclPos)[ui_idx] = (*pvclPos)[ui_idx - 1];
2070 }
2071 rvclNorm[ui_idx] = rvclNorm[ui_idx - 1];
2072 #ifdef OSG_COUNT_COMP
2073 ++g_uiSameComp;
2074 #endif
2075 }
2076 #ifdef OSG_COUNT_COMP
2077 ++g_uiTotalComp;
2078 #endif
2079 }
2080
2081 // clean up allocated memory
2082 if(i_uspan >= 0)
2083 {
2084 delete[] ppd_nu[0];
2085 delete[] ppd_nu[1];
2086 delete[] ppd_nu;
2087 }
2088 if(i_vspan >= 0)
2089 {
2090 delete[] ppd_nv[0];
2091 delete[] ppd_nv[1];
2092 delete[] ppd_nv;
2093 }
2094 delete[] apcl_temp[0];
2095 delete[] apcl_temp[1];
2096
2097 #ifdef OSG_COUNT_COMP
2098 cerr << g_uiTotalComp << " " << g_uiVSameComp << " " << g_uiSameComp << endl;
2099 #endif
2100
2101 }
2102
2103
2104 void BSplineTensorSurface::computeNormalTex(std::vector<Vec2d> & rvclUV,
2105 std::vector<Pnt3f> & rvclPos,
2106 std::vector<Vec3f> & rvclNorm,
2107 std::vector<Pnt2f> & rvclTexCoord,
2108 const std::vector<std::vector<Pnt2f> > *cpvvclTexCP)
2109 {
2110 const unsigned int cui_cnt = UInt32(rvclUV.size());
2111 unsigned int ui_idx;
2112 int i_uspan = -1;
2113 int i_vspan = -1;
2114 double **ppd_nu;
2115 double **ppd_nv;
2116 int i_k;
2117 int i_s;
2118 int i_l;
2119 DCTPVec4dmatrix vvcl_skl;
2120 DCTPVec3dmatrix vvcl_skl_eucl;
2121 Vec4d *apcl_temp[2];
2122 int i_r;
2123 double d_len;
2124 int i_u;
2125 int i_v;
2126 bool b_u_new;
2127 bool b_v_new;
2128 Vec4d *pcl_temps;
2129 double *pd_nuls;
2130 double *pd_nvkr;
2131 Vec4d *pcl_sklkl;
2132 Vec2d *pcl_temp_2d = new Vec2d[dimension_u + 1];
2133 Vec2d *pcl_temp_2dl;
2134 int i_old_u;
2135
2136 vvcl_skl.resize(2);
2137 vvcl_skl[0].resize(2);
2138 vvcl_skl[1].resize(2);
2139 vvcl_skl_eucl.resize(2);
2140 vvcl_skl_eucl[0].resize(2);
2141 vvcl_skl_eucl[1].resize(2);
2142
2143 apcl_temp[0] = new Vec4d[dimension_u + 1];
2144 apcl_temp[1] = new Vec4d[dimension_u + 1];
2145
2146 rvclPos.resize(cui_cnt);
2147 rvclNorm.resize(cui_cnt);
2148 rvclTexCoord.resize(cui_cnt);
2149
2150 for(ui_idx = 0; ui_idx < cui_cnt; ++ui_idx)
2151 {
2152 if(ui_idx == 0)
2153 {
2154 i_uspan = basis_function_u.computeDersBasisFuns(rvclUV[ui_idx][0], dimension_u, ppd_nu, 1);
2155 i_vspan = basis_function_v.computeDersBasisFuns(rvclUV[ui_idx][1], dimension_v, ppd_nv, 1);
2156 b_u_new = b_v_new = true;
2157 }
2158 else
2159 {
2160 if(osgAbs(rvclUV[ui_idx][0] - rvclUV[ui_idx - 1][0]) > DCTP_EPS)
2161 {
2162 i_old_u = i_uspan;
2163 if(i_uspan >= 0)
2164 {
2165 delete[] ppd_nu[0];
2166 delete[] ppd_nu[1];
2167 delete[] ppd_nu;
2168 }
2169 i_uspan = basis_function_u.computeDersBasisFuns(rvclUV[ui_idx][0], dimension_u, ppd_nu, 1);
2170 b_u_new = true;
2171 b_v_new = (i_old_u != i_uspan) ? true : false;
2172 }
2173 else
2174 {
2175 b_u_new = false;
2176 b_v_new = false;
2177 }
2178 if(osgAbs(rvclUV[ui_idx][1] - rvclUV[ui_idx - 1][1]) > DCTP_EPS)
2179 {
2180 if(i_vspan >= 0)
2181 {
2182 delete[] ppd_nv[0];
2183 delete[] ppd_nv[1];
2184 delete[] ppd_nv;
2185 }
2186 i_vspan = basis_function_v.computeDersBasisFuns(rvclUV[ui_idx][1], dimension_v, ppd_nv, 1);
2187 b_v_new = true;
2188 }
2189 }
2190 if( (i_uspan < 0) || (i_vspan < 0) )
2191 {
2192 rvclPos[ui_idx] = Pnt3f(0.0, 0.0, 0.0);
2193 rvclNorm[ui_idx] = Vec3f(0.0, 0.0, 0.0);
2194 }
2195 else if(b_v_new)
2196 {
2197 for(i_k = 0; i_k <= 1; ++i_k)
2198 {
2199 i_u = i_uspan - dimension_u;
2200 pcl_temps = apcl_temp[i_k];
2201
2202 for(i_s = 0; i_s <= dimension_u; ++i_s)
2203 {
2204 (*pcl_temps) = Vec4d(0.0, 0.0, 0.0, 0.0);
2205 i_v = i_vspan - dimension_v;
2206 pd_nvkr = ppd_nv[i_k];
2207
2208 for(i_r = 0; i_r <= dimension_v; ++i_r)
2209 {
2210 *pcl_temps += control_points[i_u][i_v] * *pd_nvkr;
2211 ++i_v;
2212 ++pd_nvkr;
2213 }
2214
2215 ++i_u;
2216 ++pcl_temps;
2217 }
2218
2219 pcl_sklkl = &(vvcl_skl[i_k][0]);
2220
2221 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
2222 {
2223 (*pcl_sklkl) = Vec4d(0.0, 0.0, 0.0, 0.0);
2224 pcl_temps = apcl_temp[i_k];
2225 pd_nuls = ppd_nu[i_l];
2226
2227 for(i_s = 0; i_s <= dimension_u; ++i_s)
2228 {
2229 *pcl_sklkl += *pcl_temps * *pd_nuls;
2230 ++pcl_temps;
2231 ++pd_nuls;
2232 }
2233
2234 ++pcl_sklkl;
2235 }
2236 }
2237
2238 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
2239
2240 // rvclPos[ ui_idx ] = aacl_skl[ 0 ][ 0 ];
2241 rvclPos[ui_idx][0] = vvcl_skl_eucl[0][0][0];
2242 rvclPos[ui_idx][1] = vvcl_skl_eucl[0][0][1];
2243 rvclPos[ui_idx][2] = vvcl_skl_eucl[0][0][2];
2244
2245 correctDers(rvclUV[ui_idx], vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
2246 // rvclNorm[ ui_idx ] = aacl_skl[ 0 ][ 1 ].cross( aacl_skl[ 1 ][ 0 ] );
2247 // d_len = rvclNorm[ ui_idx ].squareLength( );
2248 Vec3d tempnorm = vvcl_skl_eucl[0][1].cross(vvcl_skl_eucl[1][0]);
2249 rvclNorm[ui_idx][0] = tempnorm[0];
2250 rvclNorm[ui_idx][1] = tempnorm[1];
2251 rvclNorm[ui_idx][2] = tempnorm[2];
2252 d_len = tempnorm.squareLength();
2253
2254 if(d_len > DCTP_EPS * DCTP_EPS)
2255 {
2256 rvclNorm[ui_idx] *= 1.0 / sqrt(d_len);
2257 }
2258 else
2259 {
2260 rvclNorm[ui_idx][0] = rvclNorm[ui_idx][1] = rvclNorm[ui_idx][2] = 0.0;
2261 }
2262
2263 // calculate texture coord
2264 // rvclTexCoord[ ui_idx ].x() = rvclTexCoord[ ui_idx ].y() = 0.0;
2265 Vec2d temptexcoord(0, 0);
2266 i_u = i_uspan - dimension_u;
2267 pd_nuls = ppd_nu[0];
2268 pcl_temp_2dl = pcl_temp_2d;
2269
2270 for(i_l = 0; i_l <= dimension_u; ++i_l)
2271 {
2272 (*pcl_temp_2dl)[0] = (*pcl_temp_2dl)[1] = 0.0;
2273 i_v = i_vspan - dimension_v;
2274 pd_nvkr = ppd_nv[0];
2275
2276 for(i_k = 0; i_k <= dimension_v; ++i_k)
2277 {
2278 (*pcl_temp_2dl)[0] += (*cpvvclTexCP)[i_u][i_v].x() * *pd_nvkr;
2279 (*pcl_temp_2dl)[1] += (*cpvvclTexCP)[i_u][i_v].y() * *pd_nvkr;
2280 ++i_v;
2281 ++pd_nvkr;
2282 }
2283
2284 ++i_u;
2285 temptexcoord += *pcl_temp_2dl * *pd_nuls;
2286 ++pd_nuls;
2287 ++pcl_temp_2dl;
2288 rvclTexCoord[ui_idx][0] = temptexcoord[0];
2289 rvclTexCoord[ui_idx][1] = temptexcoord[1];
2290 }
2291 }
2292 else if(b_u_new)
2293 {
2294 for(i_k = 0; i_k <= 1; ++i_k)
2295 {
2296 pcl_sklkl = &(vvcl_skl[i_k][0]);
2297
2298 for(i_l = 0; i_l <= 1 - i_k; ++i_l)
2299 {
2300 (*pcl_sklkl)[0] = (*pcl_sklkl)[1] = (*pcl_sklkl)[2] = 0.0;
2301 pcl_temps = apcl_temp[i_k];
2302 pd_nuls = ppd_nu[i_l];
2303
2304 for(i_s = 0; i_s <= dimension_u; ++i_s)
2305 {
2306 *pcl_sklkl += *pcl_temps * *pd_nuls;
2307 ++pcl_temps;
2308 ++pd_nuls;
2309 }
2310
2311 ++pcl_sklkl;
2312 }
2313 }
2314
2315 RatSurfaceDerivs(vvcl_skl, vvcl_skl_eucl);
2316 //rvclPos[ ui_idx ] = aacl_skl[ 0 ][ 0 ];
2317 rvclPos[ui_idx][0] = vvcl_skl_eucl[0][0][0];
2318 rvclPos[ui_idx][1] = vvcl_skl_eucl[0][0][1];
2319 rvclPos[ui_idx][2] = vvcl_skl_eucl[0][0][2];
2320
2321 correctDers(rvclUV[ui_idx], vvcl_skl_eucl[0][0], vvcl_skl_eucl[1][0], vvcl_skl_eucl[0][1]);
2322 // rvclNorm[ ui_idx ] = aacl_skl[ 0 ][ 1 ].cross( aacl_skl[ 1 ][ 0 ] );
2323 // d_len = rvclNorm[ ui_idx ].squareLength( );
2324 Vec3d tempnorm = vvcl_skl_eucl[0][1].cross(vvcl_skl_eucl[1][0]);
2325 rvclNorm[ui_idx][0] = tempnorm[0];
2326 rvclNorm[ui_idx][1] = tempnorm[1];
2327 rvclNorm[ui_idx][2] = tempnorm[2];
2328 d_len = tempnorm.squareLength();
2329 if(d_len > DCTP_EPS * DCTP_EPS)
2330 {
2331 rvclNorm[ui_idx] *= 1.0 / sqrt(d_len);
2332 }
2333 else
2334 {
2335 rvclNorm[ui_idx][0] = rvclNorm[ui_idx][1] = rvclNorm[ui_idx][2] = 0.0;
2336 }
2337
2338 // calculate texture coord
2339 // rvclTexCoord[ ui_idx ][0] = rvclTexCoord[ ui_idx ][1] = 0.0;
2340 pd_nuls = ppd_nu[0];
2341 pcl_temp_2dl = pcl_temp_2d;
2342 Vec2d temptexcoord(0, 0);
2343
2344 for(i_l = 0; i_l <= dimension_u; ++i_l)
2345 {
2346 temptexcoord += *pcl_temp_2dl * *pd_nuls;
2347 ++pd_nuls;
2348 ++pcl_temp_2dl;
2349 }
2350
2351 rvclTexCoord[ui_idx][0] = temptexcoord[0];
2352 rvclTexCoord[ui_idx][1] = temptexcoord[1];
2353 #ifdef OSG_COUNT_COMP
2354 ++g_uiVSameComp;
2355 #endif
2356 }
2357 else
2358 {
2359 rvclPos[ui_idx] = rvclPos[ui_idx - 1];
2360 rvclNorm[ui_idx] = rvclNorm[ui_idx - 1];
2361 rvclTexCoord[ui_idx] = rvclTexCoord[ui_idx - 1];
2362 #ifdef OSG_COUNT_COMP
2363 ++g_uiSameComp;
2364 #endif
2365 }
2366 #ifdef OSG_COUNT_COMP
2367 ++g_uiTotalComp;
2368 #endif
2369 }
2370
2371 // clean up allocated memory
2372 if(i_uspan >= 0)
2373 {
2374 delete[] ppd_nu[0];
2375 delete[] ppd_nu[1];
2376 delete[] ppd_nu;
2377 }
2378 if(i_vspan >= 0)
2379 {
2380 delete[] ppd_nv[0];
2381 delete[] ppd_nv[1];
2382 delete[] ppd_nv;
2383 }
2384 delete[] apcl_temp[0];
2385 delete[] apcl_temp[1];
2386 delete[] pcl_temp_2d;
2387
2388 #ifdef OSG_COUNT_COMP
2389 cerr << g_uiTotalComp << " " << g_uiVSameComp << " " << g_uiSameComp << endl;
2390 #endif
2391 }
2392
2393
2394 void BSplineTensorSurface::computeTex(std::vector<Vec2d> & rvclUV,
2395 std::vector<Pnt3f> & rvclPos,
2396 std::vector<Pnt2f> & rvclTexCoord,
2397 const std::vector<std::vector<Vec2d> > *cpvvclTexCP)
2398 {
2399 const unsigned int cui_cnt = UInt32(rvclUV.size());
2400 unsigned int ui_idx;
2401 double *pd_nu;
2402 double *pd_nv;
2403 int i_span_u = -1;
2404 int i_span_v = -1;
2405 unsigned int ui_index_v;
2406 Vec4d *pcl_temp = new Vec4d[dimension_u + 1];
2407 Vec2d *pcl_ttemp = new Vec2d[dimension_u + 1];
2408 unsigned int ui_index_u;
2409 int i_l;
2410 int i_k;
2411 double *pd_nul;
2412 double *pd_nvk;
2413 bool b_u_new;
2414 bool b_v_new;
2415 Vec4d *pcl_templ;
2416 Vec2d *pcl_ttempl;
2417 int i_old_u;
2418
2419 rvclPos.resize(cui_cnt);
2420 rvclTexCoord.resize(cui_cnt);
2421
2422 for(ui_idx = 0; ui_idx < cui_cnt; ++ui_idx)
2423 {
2424 if(ui_idx == 0)
2425 {
2426 i_span_u = basis_function_u.computeAllNonzero(rvclUV[ui_idx][0], dimension_u, pd_nu);
2427 i_span_v = basis_function_v.computeAllNonzero(rvclUV[ui_idx][1], dimension_v, pd_nv);
2428 b_u_new = b_v_new = true;
2429 }
2430 else
2431 {
2432 if(osgAbs(rvclUV[ui_idx][0] - rvclUV[ui_idx - 1][0]) > DCTP_EPS)
2433 {
2434 i_old_u = i_span_u;
2435 if(i_span_u >= 0)
2436 delete[] pd_nu;
2437 i_span_u = basis_function_u.computeAllNonzero(rvclUV[ui_idx][0], dimension_u, pd_nu);
2438 b_u_new = true;
2439 b_v_new = (i_old_u != i_span_u) ? true : false;
2440 }
2441 else
2442 {
2443 b_u_new = false;
2444 b_v_new = false;
2445 }
2446 if(osgAbs(rvclUV[ui_idx][1] - rvclUV[ui_idx - 1][1]) > DCTP_EPS)
2447 {
2448 if(i_span_v >= 0)
2449 delete[] pd_nv;
2450 i_span_v = basis_function_v.computeAllNonzero(rvclUV[ui_idx][1], dimension_v, pd_nv);
2451 b_v_new = true;
2452 }
2453 }
2454
2455 if( (i_span_u < 0) || (i_span_v < 0) )
2456 {
2457 rvclPos[ui_idx] = Pnt3f(0.0, 0.0, 0.0);
2458 rvclTexCoord[ui_idx] = Pnt2f(0.0, 0.0);
2459 }
2460 else if(b_v_new)
2461 {
2462 ui_index_u = i_span_u - dimension_u;
2463 pd_nul = pd_nu;
2464 pcl_templ = pcl_temp;
2465 pcl_ttempl = pcl_ttemp;
2466 Vec4d tempposition(0.0, 0.0, 0.0, 0.0);
2467 Vec2d temptex(0, 0);
2468
2469 for(i_l = 0; i_l <= dimension_u; ++i_l)
2470 {
2471 (*pcl_templ) = Vec4d(0.0, 0.0, 0.0, 0.0);
2472 (*pcl_ttempl) = Vec2d(0.0, 0.0);
2473 ui_index_v = i_span_v - dimension_v;
2474 pd_nvk = pd_nv;
2475
2476 for(i_k = 0; i_k <= dimension_v; ++i_k)
2477 {
2478 *pcl_templ += control_points[ui_index_u][ui_index_v] * *pd_nvk;
2479 *pcl_ttempl += (*cpvvclTexCP)[ui_index_u][ui_index_v] * *pd_nvk;
2480 ++ui_index_v;
2481 ++pd_nvk;
2482 }
2483
2484 ++ui_index_u;
2485 tempposition += *pcl_templ * *pd_nul;
2486 temptex += *pcl_ttempl * *pd_nul;
2487 ++pd_nul;
2488 ++pcl_templ;
2489 ++pcl_ttempl;
2490 }
2491
2492 rvclPos[ui_idx][0] = tempposition[0] / tempposition[3];
2493 rvclPos[ui_idx][1] = tempposition[1] / tempposition[3];
2494 rvclPos[ui_idx][2] = tempposition[2] / tempposition[3];
2495 rvclTexCoord[ui_idx][0] = temptex[0];
2496 rvclTexCoord[ui_idx][1] = temptex[1];
2497 }
2498 else if(b_u_new)
2499 {
2500 ui_index_u = i_span_u - dimension_u;
2501 pd_nul = pd_nu;
2502 pcl_templ = pcl_temp;
2503 pcl_ttempl = pcl_ttemp;
2504 Vec4d tempposition(0.0, 0.0, 0.0, 0.0);
2505 Vec2d temptex(0, 0);
2506
2507 for(i_l = 0; i_l <= dimension_u; ++i_l)
2508 {
2509 ++ui_index_u;
2510 tempposition += *pcl_templ * *pd_nul;
2511 temptex += *pcl_ttempl * *pd_nul;
2512 ++pd_nul;
2513 ++pcl_templ;
2514 ++pcl_ttempl;
2515 }
2516
2517 rvclPos[ui_idx][0] = tempposition[0] / tempposition[3];
2518 rvclPos[ui_idx][1] = tempposition[1] / tempposition[3];
2519 rvclPos[ui_idx][2] = tempposition[2] / tempposition[3];
2520 rvclTexCoord[ui_idx][0] = temptex[0];
2521 rvclTexCoord[ui_idx][1] = temptex[1];
2522 #ifdef OSG_COUNT_COMP
2523 ++g_uiVSameComp;
2524 #endif
2525 }
2526 else
2527 {
2528 rvclPos[ui_idx] = rvclPos[ui_idx - 1];
2529 rvclTexCoord[ui_idx] = rvclTexCoord[ui_idx - 1];
2530 #ifdef OSG_COUNT_COMP
2531 ++g_uiSameComp;
2532 #endif
2533 }
2534 #ifdef OSG_COUNT_COMP
2535 ++g_uiTotalComp;
2536 #endif
2537 }
2538
2539 if(i_span_u >= 0)
2540 delete[] pd_nu;
2541 if(i_span_v >= 0)
2542 delete[] pd_nv;
2543 delete[] pcl_temp;
2544 delete[] pcl_ttemp;
2545
2546 #ifdef OSG_COUNT_COMP
2547 cerr << g_uiTotalComp << " " << g_uiVSameComp << " " << g_uiSameComp << endl;
2548 #endif
2549 }
2550
2551
2552 void BSplineTensorSurface::computeTexforTrimming(std::vector<Vec2d> & rvclUV,
2553 std::vector<Vec2d> & rvclTexCoord,
2554 const std::vector<std::vector<Vec2d> > *cpvvclTexCP)
2555 {
2556 const unsigned int cui_cnt = UInt32(rvclUV.size());
2557 unsigned int ui_idx;
2558 double *pd_nu;
2559 double *pd_nv;
2560 int i_span_u = -1;
2561 int i_span_v = -1;
2562 unsigned int ui_index_v;
2563 Vec2d *pcl_ttemp = new Vec2d[dimension_u + 1];
2564 unsigned int ui_index_u;
2565 int i_l;
2566 int i_k;
2567 double *pd_nul;
2568 double *pd_nvk;
2569 bool b_u_new;
2570 bool b_v_new;
2571 Vec2d *pcl_ttempl;
2572 int i_old_u;
2573
2574 rvclTexCoord.resize(cui_cnt);
2575
2576 for(ui_idx = 0; ui_idx < cui_cnt; ++ui_idx)
2577 {
2578 if(ui_idx == 0)
2579 {
2580 i_span_u = basis_function_u.computeAllNonzero(rvclUV[ui_idx][0], dimension_u, pd_nu);
2581 i_span_v = basis_function_v.computeAllNonzero(rvclUV[ui_idx][1], dimension_v, pd_nv);
2582 b_u_new = b_v_new = true;
2583 }
2584 else
2585 {
2586 if(osgAbs(rvclUV[ui_idx][0] - rvclUV[ui_idx - 1][0]) > DCTP_EPS)
2587 {
2588 i_old_u = i_span_u;
2589 if(i_span_u >= 0)
2590 delete[] pd_nu;
2591 i_span_u = basis_function_u.computeAllNonzero(rvclUV[ui_idx][0], dimension_u, pd_nu);
2592 b_u_new = true;
2593 b_v_new = (i_old_u != i_span_u) ? true : false;
2594 }
2595 else
2596 {
2597 b_u_new = false;
2598 b_v_new = false;
2599 }
2600 if(osgAbs(rvclUV[ui_idx][1] - rvclUV[ui_idx - 1][1]) > DCTP_EPS)
2601 {
2602 if(i_span_v >= 0)
2603 delete[] pd_nv;
2604 i_span_v = basis_function_v.computeAllNonzero(rvclUV[ui_idx][1], dimension_v, pd_nv);
2605 b_v_new = true;
2606 }
2607 }
2608
2609 if( (i_span_u < 0) || (i_span_v < 0) )
2610 {
2611 rvclTexCoord[ui_idx][0] = rvclTexCoord[ui_idx][1] = 0.0;
2612 }
2613 else if(b_v_new)
2614 {
2615 ui_index_u = i_span_u - dimension_u;
2616 pd_nul = pd_nu;
2617 pcl_ttempl = pcl_ttemp;
2618 rvclTexCoord[ui_idx][0] = rvclTexCoord[ui_idx][1] = 0.0;
2619
2620 for(i_l = 0; i_l <= dimension_u; ++i_l)
2621 {
2622 (*pcl_ttempl)[0] = (*pcl_ttempl)[1] = 0.0;
2623 ui_index_v = i_span_v - dimension_v;
2624 pd_nvk = pd_nv;
2625
2626 for(i_k = 0; i_k <= dimension_v; ++i_k)
2627 {
2628 *pcl_ttempl += (*cpvvclTexCP)[ui_index_u][ui_index_v] * *pd_nvk;
2629 ++ui_index_v;
2630 ++pd_nvk;
2631 }
2632
2633 ++ui_index_u;
2634 rvclTexCoord[ui_idx] += *pcl_ttempl * *pd_nul;
2635 ++pd_nul;
2636 ++pcl_ttempl;
2637 }
2638 }
2639 else if(b_u_new)
2640 {
2641 ui_index_u = i_span_u - dimension_u;
2642 pd_nul = pd_nu;
2643 pcl_ttempl = pcl_ttemp;
2644 rvclTexCoord[ui_idx][0] = rvclTexCoord[ui_idx][1] = 0.0;
2645
2646 for(i_l = 0; i_l <= dimension_u; ++i_l)
2647 {
2648 ++ui_index_u;
2649 rvclTexCoord[ui_idx] += *pcl_ttempl * *pd_nul;
2650 ++pd_nul;
2651 ++pcl_ttempl;
2652 }
2653
2654 #ifdef OSG_COUNT_COMP
2655 ++g_uiVSameComp;
2656 #endif
2657 }
2658 else
2659 {
2660 rvclTexCoord[ui_idx] = rvclTexCoord[ui_idx - 1];
2661 #ifdef OSG_COUNT_COMP
2662 ++g_uiSameComp;
2663 #endif
2664 }
2665 #ifdef OSG_COUNT_COMP
2666 ++g_uiTotalComp;
2667 #endif
2668 }
2669
2670 if(i_span_u >= 0)
2671 delete[] pd_nu;
2672 if(i_span_v >= 0)
2673 delete[] pd_nv;
2674 delete[] pcl_ttemp;
2675
2676 #ifdef OSG_COUNT_COMP
2677 cerr << g_uiTotalComp << " " << g_uiVSameComp << " " << g_uiSameComp << endl;
2678 #endif
2679 }
2680
2681
2682 void BSplineTensorSurface::correctDers(const Vec2d cclUV, const Vec3d cclPos, Vec3d &rclDU, Vec3d &rclDV)
2683 {
2684 return;
2685 // TODO: this tried to correct 0 derivs, but it was
2686 // TODO: only half-implemented plus the approach didn't really work.
2687 // TODO: implement something like this someday...
2688 #if 0
2689 if(rclDU.squareLength() < DCTP_EPS)
2690 {
2691 // du is zero
2692 Vec3d cl_temp;
2693 Vec3d cl_low;
2694 Vec3d cl_high;
2695 Vec2d cl_param;
2696 double d_step;
2697 double d_minpar;
2698 double d_maxpar;
2699 int i_err;
2700
2701 getParameterInterval_U(d_minpar, d_maxpar);
2702
2703 // step in -u direction
2704 cl_param = cclUV;
2705 cl_low = cclPos;
2706 d_step = cl_param[0] - d_minpar;
2707 while(d_step > DCTP_EPS)
2708 {
2709 d_step *= 0.5;
2710 cl_param[0] -= d_step;
2711 cl_temp = compute(cl_param, i_err);
2712 if( (cl_temp - cclPos).squareLength() > DCTP_EPS)
2713 {
2714 cl_low = cl_temp;
2715 cl_param[0] += d_step;
2716 }
2717 /* else if( ( cl_temp - cclPos ).squareLength( ) > ( cl_low - cclPos ).squareLength( ) )
2718 {
2719 cl_low = cl_temp;
2720 }*/
2721 }
2722
2723 // step in +u direction
2724 cl_param = cclUV;
2725 cl_high = cclPos;
2726 d_step = d_maxpar - cl_param[0];
2727 while(d_step > DCTP_EPS)
2728 {
2729 d_step *= 0.5;
2730 cl_param[0] += d_step;
2731 cl_temp = compute(cl_param, i_err);
2732 if( (cl_temp - cclPos).squareLength() > DCTP_EPS)
2733 {
2734 cl_high = cl_temp;
2735 cl_param[0] -= d_step;
2736 }
2737 /* else if( ( cl_temp - cclPos ).squareLength( ) > ( cl_high - cclPos ).squareLength( ) )
2738 {
2739 cl_high = cl_temp;
2740 }*/
2741 }
2742
2743 // calculate dirivative
2744 // std::cerr << "du: " << rclDU;
2745 if( (cl_high - cl_low).squareLength() > DCTP_EPS)
2746 {
2747 rclDU = cl_high - cl_low;
2748 }
2749 else
2750 {
2751 rclDU[0] = rclDU[1] = rclDU[2] = 0.0;
2752 }
2753 // std::cerr << " -> " << rclDU << std::endl;
2754 }
2755
2756 if(rclDV.squareLength() < DCTP_EPS)
2757 {
2758 // dv is zero
2759 Vec3d cl_temp;
2760 Vec3d cl_low;
2761 Vec3d cl_high;
2762 Vec2d cl_param;
2763 double d_step;
2764 double d_minpar;
2765 double d_maxpar;
2766 int i_err;
2767
2768 getParameterInterval_V(d_minpar, d_maxpar);
2769
2770 // step in -v direction
2771 cl_param = cclUV;
2772 cl_low = cclPos;
2773 d_step = cl_param[1] - d_minpar;
2774 while(d_step > DCTP_EPS)
2775 {
2776 d_step *= 0.5;
2777 cl_param[1] -= d_step;
2778 cl_temp = compute(cl_param, i_err);
2779 if( (cl_temp - cclPos).squareLength() > DCTP_EPS)
2780 {
2781 cl_low = cl_temp;
2782 cl_param[1] += d_step;
2783 }
2784 /* else if( ( cl_temp - cclPos ).squareLength( ) > ( cl_low - cclPos ).squareLength( ) )
2785 {
2786 cl_low = cl_temp;
2787 }*/
2788 }
2789
2790 // step in +v direction
2791 cl_param = cclUV;
2792 cl_high = cclPos;
2793 d_step = d_maxpar - cl_param[1];
2794 while(d_step > DCTP_EPS)
2795 {
2796 d_step *= 0.5;
2797 cl_param[1] += d_step;
2798 cl_temp = compute(cl_param, i_err);
2799 if( (cl_temp - cclPos).squareLength() > DCTP_EPS)
2800 {
2801 cl_high = cl_temp;
2802 cl_param[1] -= d_step;
2803 }
2804 /* else if( ( cl_temp - cclPos ).squareLength( ) > ( cl_high - cclPos ).squareLength( ) )
2805 {
2806 cl_high = cl_temp;
2807 }*/
2808 }
2809
2810 // calculate dirivative
2811 // std::cerr << "dv: " << rclDV;
2812 if( (cl_high - cl_low).squareLength() > DCTP_EPS)
2813 {
2814 rclDV = cl_high - cl_low;
2815 }
2816 else
2817 {
2818 rclDV[0] = rclDV[1] = rclDV[2] = 0.0;
2819 }
2820 // std::cerr << " -> " << rclDV << std::endl;
2821 }
2822 #endif
2823 }