1 /*=====================================================================
2 File: hnurbsS_sp.cpp
3 Purpose:
4 Revision: $Id: hnurbsS_sp.cpp,v 1.2 2002/05/13 21:07:46 philosophil Exp $
5 Created by: Philippe Lavoie (14 May, 1998)
6 Modified by:
7
8 Copyright notice:
9 Copyright (C) 1996-1997 Philippe Lavoie
10
11 This library is free software; you can redistribute it and/or
12 modify it under the terms of the GNU Library General Public
13 License as published by the Free Software Foundation; either
14 version 2 of the License, or (at your option) any later version.
15
16 This library is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 Library General Public License for more details.
20
21 You should have received a copy of the GNU Library General Public
22 License along with this library; if not, write to the Free
23 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 =====================================================================*/
25 #include <hnurbsS_sp.h>
26 #include <string.h>
27
28 /*!
29 */
30 namespace PLib {
31
32 /*!
33 \brief Updates the basis value in the U direction
34
35 Updates the basis value at which a control point has
36 maximal influence. It also finds where the control point
37 has maximal influence.
38
39 \warning The degrre in U of the surface must be of 3 or less.
40
41 \author Philippe Lavoie
42 \date 14 May, 1998
43 */
44 template <class T, int N>
updateMaxU()45 void HNurbsSurfaceSP<T,N>::updateMaxU() {
46 if(this->degU>3){
47 #ifdef USE_EXCEPTION
48 throw NurbsError();
49 #else
50 Error error("HNurbsSurfaceSP<T,N>::updateMaxU()") ;
51 error << "This class doesn't support surfaces having a degree of 4 and higher.\n" ;
52 error.fatal() ;
53 #endif
54 }
55 else{
56 maxU.resize(this->P.rows()) ;
57 maxAtU_.resize(this->P.rows()) ;
58 for(int i=0;i<this->P.rows();++i){
59 if(!maxInfluence(i,this->U,this->degU,maxAtU_[i]))
60 cerr << "Problem in maxInfluence U!\n" ;
61 maxU[i] = nurbsBasisFun(maxAtU_[i],i,this->degU,this->U) ;
62 }
63
64 }
65 }
66
67 /*!
68 \brief Updates the basis value in the V direction
69
70 Updates the basis value at which a control point has
71 maximal influence. It also finds where the control point
72 has maximal influence.
73
74 \warning The degrre in V of the surface must be of 3 or less.
75
76 \author Philippe Lavoie
77 \date 14 May, 1998
78 */
79 template <class T, int N>
updateMaxV()80 void HNurbsSurfaceSP<T,N>::updateMaxV() {
81 if(this->degV>3){
82 #ifdef USE_EXCEPTION
83 throw NurbsError();
84 #else
85 Error error("HNurbsSurfaceSP<T,N>::updateMaxV()") ;
86 error << "This class doesn't support surfaces having a degree of 4 and higher.\n" ;
87 error.fatal() ;
88 #endif
89 }
90 else{
91 maxV.resize(this->P.cols()) ;
92 maxAtV_.resize(this->P.cols()) ;
93 for(int i=0;i<this->P.cols();++i){
94 if(!maxInfluence(i,this->V,this->degV,maxAtV_[i]))
95 cerr << "Problem in maxInfluence V!\n" ;
96 maxV[i] = nurbsBasisFun(maxAtV_[i],i,this->degV,this->V) ;
97 }
98
99 }
100 }
101
102 /*!
103 \brief Modifies the surface point by a certain value.
104
105 \param i the row of the surface point
106 \param j the column of the surface point
107 \param a modify the surface point by this value
108
109 \warning The degree in U and V of the surface must be of 3 or less.
110
111 \author Philippe Lavoie
112 \date 14 May, 1998
113 */
114 template <class T, int N>
modSurfCPby(int i,int j,const HPoint_nD<T,N> & a)115 void HNurbsSurfaceSP<T,N>::modSurfCPby(int i, int j, const HPoint_nD<T,N>& a) {
116 this->offset(i,j) += a / (maxU[i]*maxV[j]) ;
117 if(this->baseLevel_){
118 Point_nD<T,N> vecOffset ;
119 vecOffset = this->offset(i,j).x()*this->ivec(i,j) +
120 this->offset(i,j).y()*this->jvec(i,j) +
121 this->offset(i,j).z()*this->kvec(i,j) ;
122 this->P(i,j).x() = this->baseSurf.ctrlPnts()(i,j).x()+vecOffset.x() ;
123 this->P(i,j).y() = this->baseSurf.ctrlPnts()(i,j).y()+vecOffset.y() ;
124 this->P(i,j).z() = this->baseSurf.ctrlPnts()(i,j).z()+vecOffset.z() ;
125 }
126 else
127 this->P(i,j) = this->offset(i,j) ;
128 }
129
130 /*!
131 \brief Moves the surface point only
132
133 Moves only the specified surface point. The other surface
134 points normally affected by moving this point are {\em not}
135 moved.
136
137 The point a is in the 4D homogenous space, but only
138 the x,y,z value are used. The weight is not moved by
139 this function.
140
141 \param i the row of the surface point to move
142 \param j the column of the surface point to move
143 \param a move that surface point by that amount.
144
145 \warning The degree of the curve must be of 3 or less.
146
147 \author Philippe Lavoie
148 \date 7 June, 1998
149 */
150 template <class T, int N>
modOnlySurfCPby(int i,int j,const HPoint_nD<T,N> & a)151 void HNurbsSurfaceSP<T,N>::modOnlySurfCPby(int i, int j, const HPoint_nD<T,N>& a){
152 int k ;
153
154 this->P = this->offset ;
155
156 // by definition the offset has w = 0 , but this isn't valid for
157 // the control points, increasing the w by 1, will generate a valid surface
158 if(this->baseLevel_)
159 for(k=0;k<this->P.rows();++k)
160 for(int l=0;l<this->P.cols();++l)
161 this->P(k,l).w() += T(1) ;
162
163 int sizeU, sizeV ;
164
165 sizeU = 2*this->degU+3 ;
166 if(i-this->degU-1<0) sizeU += i-this->degU-1 ;
167 if(i+this->degU+1>=this->P.rows()) sizeU -= i+this->degU+1-this->P.rows() ;
168
169 sizeV = 2*this->degV+3 ;
170 if(j-this->degV-1<0) sizeV += j-this->degV-1 ;
171 if(j+this->degV+1>=this->P.cols()) sizeV -= j+this->degV+1-this->P.cols() ;
172
173 Vector<T> u(sizeU) ;
174 Vector<T> v(sizeV) ;
175 Vector<Point_nD<T,N> > pts(sizeU*sizeV) ;
176 Vector<int> pu(sizeU*sizeV) ;
177 Vector<int> pv(sizeU*sizeV) ;
178
179 int n=0;
180 int nu = 0 ;
181 int nv = 0 ;
182 for(k=i-this->degU-1;k<=i+this->degU+1;++k){
183 if(k<0)
184 continue ;
185 if(k>=this->P.rows())
186 break ;
187 nv = 0 ;
188 for(int l=j-this->degV-1;l<=j+this->degV+1;++l){
189 if(l<0)
190 continue ;
191 if(l>=this->P.cols())
192 break ;
193 if( k == i && j==l){
194 pts[n].x() = a.x() ;
195 pts[n].y() = a.y() ;
196 pts[n].z() = a.z() ;
197 }
198 //else
199 //pts[n] = Point_nD<T,N>(0,0,0) ;
200 pu[n] = nu ;
201 pv[n] = nv ;
202 if(k==i){
203 v[nv] = maxAtV_[l] ; // only need to initialise this once
204 }
205 ++n ;
206 ++nv ;
207 }
208 u[nu] = maxAtU_[k] ;
209 ++nu ;
210 }
211
212 u.resize(nu) ;
213 v.resize(nv) ;
214 pts.resize(n) ;
215 pu.resize(n) ;
216 pv.resize(n) ;
217
218 if(NurbsSurface<T,N>::movePoint(u,v,pts,pu,pv)){
219 this->offset = this->P ;
220 // an offset shouldn't have a weight value.
221 if(this->baseLevel_)
222 for(k=0;k<this->P.rows();++k)
223 for(int l=0;l<this->P.cols();++l)
224 this->offset(k,l).w() -= T(1) ;
225 }
226 updateSurface();
227 }
228
229 /*
230 not defined for now. I'm not sure we need it/want it.
231
232 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
233 Modifies the surface point to a certain value.
234 \param i the row of the surface point
235 j the column of the surface point
236 a modify the surface point to this value
237 \return
238 \warning The degree in U and V of the surface must be of 3 or less.
239
240 \author Philippe Lavoie
241 \date 14 May, 1998
242 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
243 /*void HNurbsSurfaceSP<T,N>::modSurfCP(int i, int j, const HPoint_nD<T,N>& a) {
244 modSurfCPby(i,j,a-surfP(i,j)) ;
245 }
246 */
247
248 /*!
249 \brief Adds a level to this HNURBS surface
250
251 \param n the number of new knots between each knots.
252 \param s the multiplicity of each of these knots
253
254 \return a pointer to the new level or 0 if there was an error.
255
256 \warning returns 0, if there is already a nextlevel.
257
258 \author Philippe Lavoie
259 \date 14 May 1998
260 */
261 template <class T, int N>
addLevel(int n,int s)262 HNurbsSurfaceSP<T,N>* HNurbsSurfaceSP<T,N>::addLevel(int n, int s) {
263 HNurbsSurfaceSP<T,N> *newLevel ;
264
265 if(this->nextLevel_)
266 return 0 ;
267
268 Vector<T> newU,newV ;
269
270 this->splitUV(n,s,n,s,newU,newV) ;
271
272 newLevel = new HNurbsSurfaceSP<T,N>(this,newU,newV) ;
273
274 return newLevel ;
275 }
276
277 /*!
278 \brief Adds a level to this HNURBS surface
279
280 \param n the number of new knots between each knots.
281
282 \return a pointer to the new level or 0 if there was an error.
283 \warning returns 0, if there is already a nextlevel.
284
285 \author Philippe Lavoie
286 \date 14 May 1998
287 */
288 template <class T, int N>
addLevel()289 HNurbsSurfaceSP<T,N>* HNurbsSurfaceSP<T,N>::addLevel() {
290 HNurbsSurfaceSP<T,N> *newLevel ;
291
292 if(this->nextLevel_)
293 return 0 ;
294
295 newLevel = new HNurbsSurfaceSP<T,N>(this) ;
296
297 return newLevel ;
298 }
299
300 /*!
301 \brief Copies a HNurbs Surface and all it children
302
303 \param ns the HNurbs surface to copy
304
305 \author Philippe Lavoie
306 \date 14 May 1998
307 */
308 template <class T, int N>
copy(const HNurbsSurface<T,N> & nS)309 void HNurbsSurfaceSP<T,N>::copy(const HNurbsSurface<T,N>& nS){
310 HNurbsSurface<T,N> *levelP ;
311 levelP = nS.nextLevel() ;
312
313 NurbsSurface<T,N>::operator=(nS) ;
314 this->rU = nS.rU ;
315 this->rV = nS.rV ;
316 this->offset = nS.offset ;
317
318 updateMaxUV() ;
319
320 this->firstLevel_ = this ;
321
322 if(levelP){
323 HNurbsSurfaceSP<T,N> *newLevel ;
324 newLevel = new HNurbsSurfaceSP<T,N>(this) ;
325 newLevel->copy(*levelP) ;
326 this->nextLevel_ = newLevel ;
327 this->lastLevel_ = this->nextLevel_->lastLevel() ;
328 }
329 else{
330 this->lastLevel_ = this ;
331 }
332
333 }
334
335 /*!
336 \brief Updates the NURBS surface
337
338 Updates the NURBS surface according to the offset values
339 and its base level. You can update only one control point
340 from the surface if you specify a value for i and j or you
341 can update all the points if i0 or j0 is below 0.
342
343 \param i0 the row of the control point to update
344 \param j0 the column of the control point to update
345
346 \author Philippe Lavoie
347 \date 7 October 1997
348 */
349 template <class T, int N>
updateSurface(int i0,int j0)350 void HNurbsSurfaceSP<T,N>::updateSurface(int i0, int j0){
351 if(i0>=0 && j0>=0){
352 if(this->offset(i0,j0).x()==0.0 && this->offset(i0,j0).y()==0.0 && this->offset(i0,j0).z()==0.0)
353 return ;
354 }
355 if(this->baseLevel_){
356 if(this->initBase()){
357 this->P = this->baseSurf.ctrlPnts() ;
358 this->U = this->baseSurf.knotU() ;
359 this->V = this->baseSurf.knotV() ;
360 this->degU = this->baseSurf.degreeU() ;
361 this->degV = this->baseSurf.degreeV() ;
362 updateMaxUV() ;
363 }
364 if(i0>=0 && j0>=0){
365 Point_nD<T,N> vecOffset ;
366 vecOffset = this->offset(i0,j0).x()*this->ivec(i0,j0) +
367 this->offset(i0,j0).y()*this->jvec(i0,j0) +
368 this->offset(i0,j0).z()*this->kvec(i0,j0) ;
369 this->P(i0,j0).x() = this->baseSurf.ctrlPnts()(i0,j0).x()+vecOffset.x() ;
370 this->P(i0,j0).y() = this->baseSurf.ctrlPnts()(i0,j0).y()+vecOffset.y() ;
371 this->P(i0,j0).z() = this->baseSurf.ctrlPnts()(i0,j0).z()+vecOffset.z() ;
372 }
373 else{
374 for(int i=0;i<this->P.rows();++i)
375 for(int j=0;j<this->P.cols();++j){
376 if(this->offset(i,j).x() != 0 ||
377 this->offset(i,j).y() != 0 || this->offset(i,j).z() || 0){
378 Point_nD<T,N> vecOffset ;
379 vecOffset = this->offset(i,j).x()*this->ivec(i,j) +
380 this->offset(i,j).y()*this->jvec(i,j) +
381 this->offset(i,j).z()*this->kvec(i,j) ;
382 this->P(i,j).x() = this->baseSurf.ctrlPnts()(i,j).x()+vecOffset.x() ;
383 this->P(i,j).y() = this->baseSurf.ctrlPnts()(i,j).y()+vecOffset.y() ;
384 this->P(i,j).z() = this->baseSurf.ctrlPnts()(i,j).z()+vecOffset.z() ;
385 }
386 }
387 }
388 }
389 else{
390 if(i0>=0 && j0>=0)
391 this->P(i0,j0) = this->offset(i0,j0) ;
392 else{
393 for(int i=0;i<this->P.rows();++i)
394 for(int j=0;j<this->P.cols();++j){
395 this->P(i,j) = this->offset(i,j) ;
396 }
397 }
398 }
399
400 ++(this->updateN) ;
401 }
402
403 /*!
404 \brief Update the surface for all the levels
405
406 \param upLevel updates the levels up to this level of detail
407
408 \author Philippe Lavoie
409 \date 7 October 1997
410 */
411 template <class T, int N>
updateLevels(int upLevel)412 void HNurbsSurfaceSP<T,N>::updateLevels(int upLevel){
413 if(!okMax())
414 updateMaxUV() ;
415 if(upLevel>=0){
416 if(this->level()<=upLevel){
417 this->updateSurface() ;
418 }
419 }
420 else{
421 this->updateSurface() ;
422 }
423
424 if(upLevel>this->level() || upLevel<0){
425 if(this->nextLevel_)
426 ((HNurbsSurfaceSP<T,N>*)this->nextLevel_)->updateLevels(upLevel) ;
427 }
428 }
429
430 /*!
431 \brief Read a HNURBS surface from an input file stream.
432
433 \param fin the input file stream
434
435 \return 0 if an error occurs, 1 otherwise
436
437 \author Philippe Lavoie
438 \date 7 October 1997
439 */
440 template <class T, int N>
read(ifstream & fin)441 int HNurbsSurfaceSP<T,N>::read(ifstream &fin){
442 if(!fin) {
443 return 0 ;
444 }
445 int nu,nv,du,dv;
446 char *type ;
447 type = new char[4] ;
448 if(!fin.read(type,sizeof(char)*4)) { delete []type ; return 0 ;}
449 int r1 = strncmp(type,"hns3",4) ;
450 int r2 = strncmp(type,"hns4",4) ;
451 int r3 = strncmp(type,"hnso",4) ;
452 if(!(r1 || r2 || r3))
453 return 0 ;
454 T *p,*p2 ;
455 if(!r1 || !r2){
456 if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;}
457 if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;}
458 if(!fin.read((char*)&du,sizeof(int))) { delete []type ; return 0 ;}
459 if(!fin.read((char*)&dv,sizeof(int))) { delete []type ; return 0 ;}
460
461 this->resize(nu,nv,du,dv) ;
462
463 if(!fin.read((char*)this->U.memory(),sizeof(T)*this->U.n())) { delete []type ; return 0 ;}
464 if(!fin.read((char*)this->V.memory(),sizeof(T)*this->V.n())) { delete []type ; return 0 ;}
465
466 if(!r1){
467 p = new T[3*nu*nv] ;
468 if(!fin.read((char*)p,sizeof(T)*3*nu*nv)) { delete []type ; return 0 ;}
469 p2 = p ;
470 for(int i=0;i<nu;i++)
471 for(int j=0;j<nv;j++){
472 this->P(i,j).x() = *(p++) ;
473 this->P(i,j).y() = *(p++) ;
474 this->P(i,j).z() = *(p++) ;
475 this->P(i,j).w() = 1.0 ;
476 }
477 delete []p2 ;
478 }
479 else{
480 p = new T[4*nu*nv] ;
481 if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
482 p2 = p ;
483 for(int i=0;i<nu;i++)
484 for(int j=0;j<nv;j++){
485 this->P(i,j).x() = *(p++) ;
486 this->P(i,j).y() = *(p++) ;
487 this->P(i,j).z() = *(p++) ;
488 this->P(i,j).w() = *(p++) ;
489 }
490 delete []p2 ;
491 }
492 this->offset = this->P ;
493 this->updateSurface() ;
494 }
495 else { // reading the offset information
496 int ru,rv ;
497 if(!fin.read((char*)&ru,sizeof(int))) { delete []type ; return 0 ;}
498 if(!fin.read((char*)&rv,sizeof(int))) { delete []type ; return 0 ;}
499 this->rU.resize(ru) ;
500 this->rV.resize(rv) ;
501 if(this->rU.n()>0)
502 if(!fin.read((char*)this->rU.memory(),sizeof(T)*this->rU.n())) { delete []type ; return 0 ;}
503 if(this->rV.n()>0)
504 if(!fin.read((char*)this->rV.memory(),sizeof(T)*this->rV.n())) { delete []type ; return 0 ;}
505
506 if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;}
507 if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;}
508
509 p = new T[4*nu*nv] ;
510 if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
511 p2 = p ;
512 this->offset.resize(nu,nv) ;
513 for(int i=0;i<nu;i++)
514 for(int j=0;j<nv;j++){
515 this->offset(i,j).x() = *(p++) ;
516 this->offset(i,j).y() = *(p++) ;
517 this->offset(i,j).z() = *(p++) ;
518 this->offset(i,j).w() = *(p++) ;
519 }
520 delete []p2 ;
521 --(this->baseUpdateN) ;
522 this->updateSurface() ;
523 }
524
525 // now we must read all the level information (if any)
526 char *ptxt ;
527 ptxt = new char[7] ;
528 // the following line is used so that purify doesn't give a warning
529 ptxt[0] = ptxt[1] = ptxt[2] = ptxt[3] = ptxt[4] = ptxt[5] = ptxt[6] = 0 ;
530 int mark = fin.tellg() ;
531 if(!fin.read(ptxt,sizeof(char)*5))
532 { delete []ptxt ; delete []type ; return 1 ; }
533 char *pat ;
534 pat = strstr(ptxt,"level") ;
535 if(pat){
536 // insert a new level
537 HNurbsSurfaceSP<T,N> *newPatch = new HNurbsSurfaceSP<T,N>(this) ;
538 if(newPatch){
539 if(!newPatch->read(fin))
540 return 0 ;
541 }
542 else
543 return 0 ;
544 }
545 else{
546 // something else is after the HNurbs
547 // restore were the file was before the read
548 fin.seekg(mark);
549 }
550
551 delete []ptxt ;
552 delete []type ;
553 return 1 ;
554 }
555
556 } // end namespace
557