00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef PLIB_NURBS_HNURBS_SOURCE
00027 #define PLIB_NURBS_HNURBS_SOURCE
00028
00029 #include <hnurbs.h>
00030
00031 HNurbsCurveNode::HNurbsCurveNode():u0(u0_),u1(u1_){
00032 prev = 0 ;
00033 next = 0 ;
00034 curve = new NurbsCurve ;
00035 u0 = 0 ;
00036 u1 = 1 ;
00037 uD = 1 ;
00038 }
00039
00040 HNurbsCurveNode::HNurbsCurveNode(const NurbsCurve &c, T uS, T uE):u0(u0_),u1(u1_){
00041 prev = 0 ;
00042 next = 0 ;
00043 curve = new NurbsCurve(c) ;
00044 u0 = uS ;
00045 u1 = uE ;
00046 uD = uE-uS ;
00047 }
00048
00049 HPoint_nD<T,N> HNurbsCurveNode::operator()(T u) const {
00050 if(u<u0 || u>u1)
00051 return HPoint_nD<T,N>(0,0,0,0) ;
00052
00053 u -= u0_ ;
00054 u /= uD ;
00055 return (*curve)(u) ;
00056 }
00057
00058 void HNurbsCurveNode::deriveAt(T u, int d, Vector< Point_nD<T,N> >& ders) const {
00059 ders.resize(d+1) ;
00060 ders.reset(0) ;
00061 if(u<u0 || u>u1)
00062 return ;
00063
00064 u -= u0_ ;
00065 u /= uD ;
00066 curve->deriveAt(u,d,ders) ;
00067 }
00068
00069 void HNurbsCurveNode::deriveAt(T u, int d, Vector< HPoint_nD<T,N> >& ders) const {
00070 ders.resize(d+1) ;
00071 ders.reset(0) ;
00072 if(u<u0 || u>u1)
00073 return ;
00074
00075 u -= u0_ ;
00076 u /= uD ;
00077 curve->deriveAt(u,d,ders) ;
00078 }
00079
00080
00081 HNurbsCurve::HNurbsCurve(){
00082 first = 0 ;
00083 last = 0 ;
00084 }
00085
00086 void HNurbsCurve::add(const NurbsCurve& curve, T uS, T uE) {
00087 HNurbsCurveNode *nC ;
00088 nC = new HNurbsCurveNode(curve,uS,uE) ;
00089 nC->prev = last ;
00090 if(last)
00091 last->next = nC ;
00092 if(!first)
00093 first = nC ;
00094 last = nC ;
00095 }
00096
00097 void HNurbsCurve::remove(void){
00098 HNurbsCurveNode *tC ;
00099 if(last){
00100 tC = last ;
00101 last = last->prev;
00102 last->next = 0 ;
00103 delete tC ;
00104 }
00105 }
00106
00107 HPoint_nD<T,N> HNurbsCurve::operator()(T u) const{
00108 HNurbsCurveNode *c ;
00109 HPoint_nD<T,N> result(0,0,0,0) ;
00110 HPoint_nD<T,N> temp ;
00111 Point_nD<T,N> temp2 ;
00112 c = first ;
00113 while(c){
00114 temp = (*c)(u) ;
00115 if(temp.w()!= 0){
00116 temp2 = project(temp) ;
00117 result.x() += temp2.x() ;
00118 result.y() += temp2.y() ;
00119 result.z() += temp2.z() ;
00120 }
00121 c = c->next ;
00122 }
00123 result.w() = 1.0 ;
00124 return result ;
00125 }
00126
00127 void HNurbsCurve::deriveAt(T u,int d, Vector< HPoint_nD<T,N> >& ders) const{
00128 HNurbsCurveNode *c ;
00129 Vector< HPoint_nD<T,N> > dTemp(d+1) ;
00130
00131 ders.resize(d+1) ;
00132 ders.reset(0) ;
00133 c = first ;
00134 while(c){
00135 c->deriveAt(u,d,dTemp) ;
00136 ders += dTemp ;
00137 c = c->next ;
00138 }
00139 }
00140 void HNurbsCurve::deriveAt(T u,int d, Vector< Point_nD<T,N> >& ders) const{
00141 HNurbsCurveNode *c ;
00142 Vector< Point_nD<T,N> > dTemp(d+1) ;
00143
00144 ders.resize(d+1) ;
00145 ders.reset(0) ;
00146 c = first ;
00147 while(c){
00148 c->deriveAt(u,d,dTemp) ;
00149 ders += dTemp ;
00150 c = c->next ;
00151 }
00152 }
00153
00154
00155
00156 void HNurbsCurve::draw(Image_Color& img, const Color& col) const {
00157 const int n=100 ;
00158 T dU = 1.0/T(n) ;
00159 T uP = 0 ;
00160 HPoint_nD<T,N> p,c ;
00161 Point_nD<T,N> p3,c3 ;
00162
00163 p = point4D(uP) ;
00164 for(T u=dU; u<1.0; u += dU){
00165 c = point4D(u) ;
00166 p3 = project(p) ;
00167 c3 = project(c) ;
00168 img.drawLine((int)p3.y(),(int)p3.x(),(int)c3.y(),(int)c3.x(),col) ;
00169 p = c ;
00170 }
00171 c = point4D(1.0) ;
00172 p3 = project(p) ;
00173 c3 = project(c) ;
00174 img.drawLine((int)p3.y(),(int)p3.x(),(int)c3.y(),(int)c3.x(),col) ;
00175 }
00176
00177 void HNurbsCurve::draw(Image_UBYTE &img, unsigned char col) const {
00178 const int n=100 ;
00179 T dU = 1.0/T(n) ;
00180 HPoint_nD<T,N> p,c ;
00181 Point_nD<T,N> p3,c3 ;
00182
00183 p3 = point3D(0.0) ;
00184 for(T u=dU;u<1.0;u+=dU){
00185 c3 = point3D(u) ;
00186 img.drawLine((int)p3.y(),(int)p3.x(),(int)c3.y(),(int)c3.x(),col) ;
00187 p3 = c3 ;
00188 }
00189 c3 = point3D(1.0);
00190 img.drawLine((int)p3.y(),(int)p3.x(),(int)c3.y(),(int)c3.x(),col) ;
00191 }
00192
00193 const int MaxRandom = 32768 ;
00194
00195
00196 void HNurbsCurve::reset(){
00197 HNurbsCurveNode* t ;
00198 while(last){
00199 t = last->prev ;
00200 if(last->curve)
00201 delete last->curve ;
00202 delete last ;
00203 last = t ;
00204 }
00205 }
00206
00207
00208
00209 void setRandomVector(Vector_INT& rV, int rangeSize, int rangeOffset, int startAt=0, int finishAt=-1){
00210 if(finishAt<0)
00211 finishAt= rV.n-1 ;
00212 if(rV.n<=finishAt){
00213 cerr << "You need to set the random vector's size to the proper value PRIOR to calling this function!\n" ;
00214 while(1){;}
00215 }
00216
00217
00218 int i ;
00219 for(i=startAt;i<=finishAt;++i){
00220 double r = double(rand()) ;
00221 int l ;
00222 r /= double(MaxRandom) ;
00223 r *= double(rangeSize) ;
00224 --rangeSize ;
00225 l = int(r) ;
00226 for(int j=startAt;j<i;++j){
00227 if(rV[j]<=l){
00228 ++l ;
00229 }
00230 else{
00231 int k=j ;
00232 for(j=i;j>k;--j){
00233 rV[j] = rV[j-1] ;
00234 }
00235 rV[k] = l ;
00236 break ;
00237 }
00238 }
00239 }
00240
00241 for(i=0;i<=finishAt-startAt;i++)
00242 rV[i] += rangeOffset ;
00243 }
00244
00245 void HNurbsCurve::interpolate(const Vector< Point_nD<T,N> > &Pts, int deg, T acceptError, int nSample, int maxTries, int nInitPoints, int nPoints){
00246
00247 if(nInitPoints<0) nInitPoints = deg*2 ;
00248 if(nPoints<0) nPoints = deg*2 ;
00249
00250 if(nPoints<deg+1 || nInitPoints < deg+1){
00251 cerr << "Using a number of points smaller than deg+1!\n" ;
00252 while(1){;}
00253 }
00254 if(Pts.n<deg+1){
00255 cerr << "Not enough points to interpolate !\n" ;
00256 while(1){;}
00257 }
00258
00259 int i,j,k ;
00260
00261 NurbsCurve sampleC,minC ;
00262 Vector< Point_nD<T,N> > sampleP ;
00263 Vector_INT sampleI ;
00264 Vector<T> error ;
00265 T minError = -1.0 ;
00266
00267 sampleP.resize(nInitPoints) ;
00268 sampleI.resize(nInitPoints) ;
00269 error.resize(Pts.n) ;
00270 srand(123) ;
00271
00272 for(k=0;k<nSample;++k){
00273 sampleI[0] = 0 ;
00274 sampleI[sampleP.n-1] = Pts.n-1 ;
00275 for(i=1;i<sampleP.n-1;++i){
00276 int ok = 0 ;
00277 while(!ok){
00278 ok = 1 ;
00279 double r = double(rand()) ;
00280 r /= double(MaxRandom) ;
00281 r *= Pts.n-2 ;
00282 r += 1 ;
00283 sampleI[i] = int(r) ;
00284 for(j=0;j<i;++j){
00285 if(sampleI[j] == sampleI[i]){
00286 ok = 0 ;
00287 break ;
00288 }
00289 }
00290 }
00291 }
00292 sampleI.qSort() ;
00293 for(i=0;i<sampleP.n;++i)
00294 sampleP[i] = Pts[sampleI[i]] ;
00295 sampleC.globalInterp(sampleP,deg) ;
00296 T u = 0 ;
00297 for(i=0;i<Pts.n;++i){
00298 T e = sampleC.minDist2(Pts[i],u) ;
00299 error[i] = e ;
00300 }
00301 error.qSort() ;
00302 if(minError>error[error.n/2] || minError<0){
00303 minError = error[error.n/2] ;
00304 minC = sampleC ;
00305 }
00306 }
00307 reset() ;
00308 add(minC,0,1) ;
00309
00310
00311
00312
00313 T maxError ;
00314 Vector_INT regionI(Pts.n) ;
00315 Vector<T> regionU(Pts.n) ;
00316 sampleP.resize(nPoints) ;
00317 sampleI.resize(nPoints) ;
00318 Vector<T> errorM,errorT ;
00319 int tryN = 0;
00320 while(tryN < maxTries){
00321
00322 T minU,maxU,tempError,u ;
00323 int maxI = 0 ;
00324
00325 int index = 0 ;
00326 regionI[0] = 0 ;
00327 regionU[0] = 0.0 ;
00328 maxError = -1.0 ;
00329 u = 0 ;
00330
00331 for(i=1;i<Pts.n-1;++i){
00332 tempError = minDist2(Pts[i],u) ;
00333 if(tempError<acceptError){
00334 ++index ;
00335 regionI[index] = i ;
00336 regionU[index] = u ;
00337 continue ;
00338 }
00339 if(tempError>maxError){
00340 maxError = tempError ;
00341 maxI = index ;
00342 }
00343 }
00344 ++index ;
00345 regionI[index] = Pts.n-1 ;
00346 regionU[index] = 1.0 ;
00347 ++index ;
00348 if(maxError<acceptError){
00349 tryN = maxTries ;
00350 break ;
00351 }
00352
00353
00354
00355 int regionS = maxI ;
00356 int regionE = maxI+1 ;
00357
00358 while((regionI[regionE]-regionI[regionS]) < 2*nPoints){
00359 if(regionS>0) --regionS ;
00360 if(regionE<index-1) ++regionE ;
00361 }
00362
00363 int regionSize = regionI[regionE]-regionI[regionS] ;
00364 int regionStart = regionI[regionS] ;
00365
00366 error.resize(regionSize) ;
00367 errorM.resize(regionSize) ;
00368
00369 minError = -1.0 ;
00370 for(k=0;k<nSample;++k){
00371 sampleI[0] = regionStart ;
00372 sampleI[sampleP.n-1] = regionI[regionE] ;
00373 if(regionSize==nPoints){
00374 k=nSample ;
00375 for(i=1;i<sampleP.n-1;++i)
00376 sampleI[i] = regionStart+i ;
00377 }
00378 else{
00379 for(i=1;i<sampleP.n-1;++i){
00380 int ok = 0 ;
00381 while(!ok){
00382 ok = 1 ;
00383 double r = double(rand()) ;
00384 r /= double(MaxRandom) ;
00385 r *= regionSize-2 ;
00386 r += 1 ;
00387 sampleI[i] = int(r)+regionStart ;
00388 for(j=0;j<i;++j){
00389 if(sampleI[j] == sampleI[i]){
00390 ok = 0 ;
00391 break ;
00392 }
00393 }
00394 }
00395 }
00396 }
00397 sampleI.qSort() ;
00398 u = regionU[regionS] ;
00399 for(i=1;i<sampleP.n-1;++i){
00400
00401 sampleP[i] = Pts[sampleI[i]] - point3D(u) ;
00402 }
00403 sampleP[0] = Point_nD<T,N>(0,0,0) ;
00404 sampleP[sampleP.n-1] = Point_nD<T,N>(0,0,0) ;
00405 sampleC.globalInterp(sampleP,deg) ;
00406 u = regionU[regionS] ;
00407 for(i=0;i<error.n;++i){
00408
00409 Point_nD<T,N> closest = point3D(u) ;
00410 T u2 = u ;
00411 u2 -= regionU[regionS] ;
00412 u2 /= regionU[regionE]-regionU[regionS] ;
00413 if(u2<0 || u2> 1.0) {
00414 cerr << "Error in setting up u2!\n" ;
00415 while(1) {; }
00416 }
00417 error[i] = abs2(closest+project(sampleC(u2))-Pts[regionStart+i]) ;
00418 }
00419 error.qSort() ;
00420 if(minError>error[error.n/2] || minError<0){
00421 minError = error[error.n/2] ;
00422 minC = sampleC ;
00423 minU = regionU[regionS] ;
00424 maxU = regionU[regionE] ;
00425 errorM = error ;
00426 }
00427 }
00428 add(minC,regionU[regionS],regionU[regionE]) ;
00429
00430 errorT.resize(errorM.n) ;
00431 u = regionU[regionS] ;
00432 for(i=0;i<errorT.n;++i) {
00433 T closestE = minDist2(Pts[i+regionStart],u) ;
00434 errorT[i] = closestE ;
00435
00436 }
00437 errorT.qSort() ;
00438
00439
00440
00441 ++tryN ;
00442 }
00443
00444 }
00445
00446
00447 #endif // #define PLIB_NURBS_HNURBS_SOURCE