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
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "G4ProjectedSurface.hh"
00037
00038 G4int G4ProjectedSurface::Splits=0;
00039
00040 G4ProjectedSurface::G4ProjectedSurface()
00041 : G4Surface(), ctl_points(0), dir(0), u_knots(0), v_knots(0),
00042 projected_list(0), bezier_list(0), new_knots(0), ord(0),
00043 lower(0), upper(0), oslo_m(0)
00044 {
00045 distance = 0;
00046 order[0] = order[1] = 0;
00047 }
00048
00049
00050 G4ProjectedSurface::~G4ProjectedSurface()
00051 {
00052 delete u_knots;
00053 delete v_knots;
00054 delete ctl_points;
00055
00056 G4OsloMatrix* temp_oslo;
00057 if(oslo_m!=(G4OsloMatrix*)0)
00058 {
00059 while(oslo_m->GetNextNode() != oslo_m)
00060 {
00061 temp_oslo = oslo_m;
00062 oslo_m = oslo_m->GetNextNode();
00063
00064 delete temp_oslo;
00065 }
00066
00067 delete oslo_m;
00068 }
00069
00070 delete bbox;
00071 }
00072
00073 G4ProjectedSurface::G4ProjectedSurface(const G4ProjectedSurface&)
00074 : G4Surface(), ctl_points(0), dir(0), u_knots(0), v_knots(0),
00075 projected_list(0), bezier_list(0), new_knots(0), ord(0),
00076 lower(0), upper(0), oslo_m(0)
00077 {
00078 order[0] = 0;
00079 order[1] = 0;
00080 }
00081
00082 void G4ProjectedSurface::CopySurface()
00083
00084
00085 {
00086 G4BezierSurface *bez = new G4BezierSurface();
00087
00088 bez->SetDistance(distance);
00089 bez->PutOrder(0, order[0]);
00090 bez->PutOrder(1, order[1]);
00091 bez->Dir(dir);
00092 bez->u_knots = new G4KnotVector(*u_knots);
00093 bez->v_knots = new G4KnotVector(*v_knots);
00094 bez->ctl_points = new G4ControlPoints(*ctl_points);
00095
00096 bezier_list->AddSurface(bez);
00097 }
00098
00099
00100 void G4ProjectedSurface::CalcBBox()
00101 {
00102
00103
00104
00105
00106
00107
00108
00109 G4double box_minx,box_miny,box_maxx,box_maxy;
00110 box_minx = kInfinity;
00111 box_miny = kInfinity;
00112 box_maxx = -kInfinity;
00113 box_maxy = -kInfinity;
00114
00115 G4double bminx,bminy,bmaxx,bmaxy,tmpx,tmpy;
00116 bminx = box_minx; bminy = box_miny;
00117 bmaxx = box_maxx; bmaxy = box_maxy;
00118
00119 for(register G4int a = ctl_points->GetRows()-1; a>=0;a--)
00120 for(register G4int b = ctl_points->GetCols()-1; b>=0;b--)
00121 {
00122
00123
00124
00125 G4Point3D tmp = ctl_points->Get3D(a,b);
00126
00127 tmpx = tmp.x(); tmpy = tmp.y();
00128 if(bminx > tmpx) box_minx=tmpx;
00129 if(bmaxx < tmpx) box_maxx=tmpx;
00130 if(bminy > tmpy) box_miny=tmpy;
00131 if(bmaxy < tmpy) box_maxy=tmpy;
00132 }
00133
00134 G4Point3D box_min(box_minx,box_miny,0.);
00135 G4Point3D box_max(box_maxx,box_maxy,0.);
00136
00137 delete bbox;
00138 bbox = new G4BoundingBox3D(box_min, box_max);
00139 }
00140
00141
00142 void G4ProjectedSurface::ConvertToBezier(G4SurfaceList& proj_list,
00143 G4SurfaceList& bez_list)
00144 {
00145 projected_list = &proj_list;
00146 bezier_list = &bez_list;
00147
00148
00149
00150 if(CheckBezier())
00151 {
00152
00153
00154 CopySurface();
00155
00156
00157
00158 G4BezierSurface* bez_ptr = (G4BezierSurface*)bezier_list->GetLastSurface();
00159
00160
00161 bez_ptr->ClipSurface();
00162 G4double dMin = bez_ptr->SMin();
00163 G4double dMax = bez_ptr->SMax();
00164 G4double dMaxMinusdMin = dMax - dMin;
00165
00166 if(( dMaxMinusdMin > kCarTolerance ))
00167 {
00168 if( dMaxMinusdMin > 0.8 )
00169 {
00170
00171
00172
00173
00174
00175
00176
00177
00178 dir = bez_ptr->dir;
00179 bezier_list->RemoveSurface(bez_ptr);
00180
00181 SplitNURBSurface();
00182 return;
00183
00184 }
00185 else
00186 if( dMin > 0.0 || dMax < 0.0 )
00187 {
00188
00189
00190
00191 bezier_list->RemoveSurface(bez_ptr);
00192 return;
00193 }
00194 }
00195 else
00196 if(dMaxMinusdMin < kCarTolerance && dMaxMinusdMin > -kCarTolerance)
00197 {
00198 bezier_list->RemoveSurface(bez_ptr);
00199 return;
00200 }
00201
00202 bez_ptr->LocalizeClipValues();
00203 bez_ptr->SetValues();
00204
00205
00206 bez_ptr->ChangeDir();
00207 bez_ptr->ClipSurface();
00208
00209
00210
00211 dMin = bez_ptr->SMin();
00212 dMax = bez_ptr->SMax();
00213 dMaxMinusdMin = dMax-dMin;
00214
00215 if((dMaxMinusdMin > kCarTolerance ))
00216
00217 {
00218 if( (dMaxMinusdMin) > 0.8 )
00219 {
00220
00221 dir = bez_ptr->dir;
00222
00223
00224 bezier_list->RemoveSurface(bez_ptr);
00225 SplitNURBSurface();
00226 return;
00227
00228 }
00229 else
00230 if( dMin > 1.0 || dMax < 0.0 )
00231 {
00232
00233 bezier_list->RemoveSurface(bez_ptr);
00234 return;
00235 }
00236 }
00237 else
00238 if(dMaxMinusdMin < kCarTolerance && dMaxMinusdMin > -kCarTolerance)
00239 {
00240 bezier_list->RemoveSurface(bez_ptr);
00241 return;
00242 }
00243
00244 bez_ptr->LocalizeClipValues();
00245 bez_ptr->SetValues();
00246 bez_ptr->CalcAverage();
00247 }
00248 else
00249 {
00250
00251
00252
00253 SplitNURBSurface();
00254 }
00255 }
00256
00257
00258 G4int G4ProjectedSurface::CheckBezier()
00259 {
00260
00261
00262
00263
00264
00265
00266 if( u_knots->GetSize() > (2.0 * GetOrder(ROW)))
00267 {dir=0;return 0;}
00268
00269 if( v_knots->GetSize() > (2.0 * GetOrder(COL)))
00270 {dir=1;return 0;}
00271
00272 return 1;
00273 }
00274
00275
00276 void G4ProjectedSurface::SplitNURBSurface()
00277 {
00278
00279
00280
00281
00282
00283 register G4double value;
00284 register G4int i;
00285 register G4int k_index=0;
00286 register G4ProjectedSurface *srf1, *srf2;
00287 register G4int nr,nc;
00288
00289 if ( dir == ROW )
00290 {
00291 value = u_knots->GetKnot((u_knots->GetSize()-1)/2);
00292
00293 for( i = 0; i < u_knots->GetSize(); i++)
00294 if( (std::abs(value - u_knots->GetKnot(i))) < kCarTolerance )
00295 {
00296 k_index = i;
00297 break;
00298 }
00299
00300 if ( k_index == 0)
00301 {
00302 value = ( value + u_knots->GetKnot(u_knots->GetSize() -1))/2.0;
00303 k_index = GetOrder(ROW);
00304 }
00305
00306 new_knots = u_knots->MultiplyKnotVector(GetOrder(ROW), value);
00307
00308 ord = GetOrder(ROW);
00309 CalcOsloMatrix();
00310
00311 srf1 = new G4ProjectedSurface(*this);
00312
00313
00314 srf1->dir=COL;
00315
00316 new_knots->ExtractKnotVector(srf1->u_knots,
00317 k_index + srf1->GetOrder(ROW),0);
00318
00319 nr= srf1->v_knots->GetSize() - srf1->GetOrder(COL);
00320 nc= srf1->u_knots->GetSize() - srf1->GetOrder(ROW);
00321
00322 delete srf1->ctl_points;
00323 srf1->ctl_points= new G4ControlPoints(2, nr, nc);
00324
00325 srf2 = new G4ProjectedSurface(*this);
00326
00327
00328 srf2->dir = COL;
00329
00330 new_knots->ExtractKnotVector(srf2->u_knots,
00331 new_knots->GetSize(), k_index);
00332
00333 nr= srf2->v_knots->GetSize() - srf2->GetOrder(COL);
00334 nc= srf2->u_knots->GetSize() - srf2->GetOrder(ROW);
00335
00336 delete srf2->ctl_points;
00337 srf2->ctl_points = new G4ControlPoints(2, nr, nc);
00338
00339 lower = 0;
00340 upper = k_index;
00341 MapSurface(srf1);
00342
00343 lower = k_index;
00344 upper = new_knots->GetSize() - srf2->GetOrder(ROW);
00345 MapSurface(srf2);
00346 }
00347 else
00348 {
00349 value = v_knots->GetKnot((v_knots->GetSize() -1)/2);
00350
00351 for( i = 0; i < v_knots->GetSize(); i++)
00352 if( (std::abs(value - v_knots->GetKnot(i))) < kCarTolerance )
00353 {
00354 k_index = i;
00355 break;
00356 }
00357
00358 if ( k_index == 0)
00359 {
00360 value = ( value + v_knots->GetKnot(v_knots->GetSize() -1))/2.0;
00361 k_index = GetOrder(COL);
00362 }
00363
00364 new_knots = v_knots->MultiplyKnotVector( GetOrder(COL), value );
00365 ord = GetOrder(COL);
00366
00367 CalcOsloMatrix();
00368
00369 srf1 = new G4ProjectedSurface(*this);
00370
00371
00372 srf1->dir = ROW;
00373
00374 new_knots->ExtractKnotVector(srf1->v_knots,
00375 k_index + srf1->GetOrder(COL), 0);
00376
00377 nr = srf1->v_knots->GetSize() - srf1->GetOrder(COL);
00378 nc = srf1->u_knots->GetSize() - srf1->GetOrder(ROW);
00379
00380 delete srf1->ctl_points;
00381 srf1->ctl_points = new G4ControlPoints(2, nr, nc);
00382
00383 srf2 = new G4ProjectedSurface(*this);
00384
00385
00386 srf2->dir = ROW;
00387
00388 new_knots->ExtractKnotVector(srf2->v_knots, new_knots->GetSize(), k_index);
00389
00390 nr = srf2->v_knots->GetSize() - srf2->GetOrder(COL);
00391 nc = srf2->u_knots->GetSize() - srf2->GetOrder(ROW);
00392
00393 delete srf2->ctl_points;
00394 srf2->ctl_points = new G4ControlPoints(2,nr, nc);
00395
00396 lower = 0;
00397 upper = k_index;
00398 MapSurface(srf1);
00399
00400 lower = k_index;
00401 upper = new_knots->GetSize() - srf2->GetOrder(COL);
00402 MapSurface(srf2);
00403 }
00404
00405
00406 G4int col_size = srf1->ctl_points->GetCols();
00407 G4int row_size = srf1->ctl_points->GetRows();
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421 G4Point3D pt1 = srf1->ctl_points->Get3D(0,0);
00422 G4Point3D pt2 = srf1->ctl_points->Get3D(0,col_size-1);
00423 G4Point3D pt3 = srf1->ctl_points->Get3D(row_size-1,0);
00424
00425
00426 G4double pointDist1 = pt1.distance2(pt2);
00427 G4double pointDist2 = pt1.distance2(pt3);
00428
00429
00430
00431 if(pointDist1 > kCarTolerance && pointDist2 > kCarTolerance)
00432 projected_list->AddSurface(srf1);
00433 else
00434 delete srf1;
00435
00436 col_size = srf2->ctl_points->GetCols();
00437 row_size = srf2->ctl_points->GetRows();
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451 pt1 = srf2->ctl_points->Get3D(0,0);
00452 pt2 = srf2->ctl_points->Get3D(0,col_size-1);
00453 pt3 = srf2->ctl_points->Get3D(row_size-1,0);
00454
00455
00456 pointDist1 = pt1.distance2(pt2);
00457 pointDist2 = pt1.distance2(pt3);
00458
00459
00460 if(pointDist1 > kCarTolerance && pointDist2 > kCarTolerance)
00461 projected_list->AddSurface(srf2);
00462 else
00463 delete srf2;
00464
00465 delete new_knots;
00466
00467 Splits++;
00468 }
00469
00470 void G4ProjectedSurface::CalcOsloMatrix()
00471 {
00472
00473
00474
00475
00476
00477 register G4KnotVector *ah;
00478 static G4KnotVector *newknots;
00479 register G4int i;
00480 register G4int j;
00481 register G4int mu, muprim;
00482 register G4int v, p;
00483 register G4int iu, il, ih, n1;
00484 register G4int ahi;
00485 register G4double beta1;
00486 register G4double tj;
00487
00488 ah = new G4KnotVector(ord*(ord + 1)/2);
00489
00490 newknots = new G4KnotVector(ord * 2 );
00491
00492 n1 = new_knots->GetSize() - ord;
00493 mu = 0;
00494
00495 if(oslo_m!=(G4OsloMatrix*)0)
00496 {
00497 G4OsloMatrix* tmp;
00498 while(oslo_m!=oslo_m->GetNextNode())
00499 {
00500 tmp=oslo_m->GetNextNode();
00501 delete oslo_m;
00502 oslo_m=tmp;
00503 }
00504 }
00505
00506 delete oslo_m;
00507
00508 oslo_m = new G4OsloMatrix();
00509 register G4OsloMatrix* o_ptr = oslo_m;
00510
00511 register G4KnotVector* old_knots;
00512
00513 if(dir)
00514 old_knots = v_knots;
00515 else
00516 old_knots = u_knots;
00517
00518 for (j = 0; j < n1; j++)
00519 {
00520 if ( j != 0 )
00521 {
00522 oslo_m->SetNextNode(new G4OsloMatrix());
00523 oslo_m = oslo_m->GetNextNode();
00524 }
00525
00526
00527 while ( (new_knots->GetKnot(j) - old_knots->GetKnot(mu + 1)) >
00528 kCarTolerance )
00529 mu = mu + 1;
00530
00531 i = j + 1;
00532 muprim = mu;
00533
00534 while ( ((std::abs(new_knots->GetKnot(i) - old_knots->GetKnot(muprim))) <
00535 kCarTolerance) && i < (j + ord) )
00536 {
00537 i++;
00538 muprim--;
00539 }
00540
00541 ih = muprim + 1;
00542
00543 for (v = 0, p = 1; p < ord; p++)
00544 {
00545
00546 if ( (std::abs((new_knots->GetKnot(j + p)) - (old_knots->GetKnot(ih)))) <
00547 kCarTolerance )
00548 ih++;
00549 else
00550 newknots->PutKnot(++v - 1,new_knots->GetKnot(j + p));
00551 }
00552
00553 ahi = AhIndex(0, ord - 1,ord);
00554 ah->PutKnot(ahi, 1.0);
00555
00556 for (p = 1; p <= v; p++)
00557 {
00558 beta1 = 0.0;
00559 tj = newknots->GetKnot(p-1);
00560
00561 if (p - 1 >= muprim)
00562 {
00563 beta1 = AhIndex(p - 1, ord - muprim,ord);
00564 beta1 = ((tj - old_knots->GetKnot(0)) * beta1) /
00565 (old_knots->GetKnot(p + ord - v) - old_knots->GetKnot(0));
00566 }
00567
00568 i = muprim - p + 1;
00569 il = Amax (1, i);
00570 i = n1 - 1 + v - p;
00571 iu = Amin (muprim, i);
00572
00573 for (i = il; i <= iu; i++)
00574 {
00575 register G4double d1, d2;
00576 register G4double beta;
00577
00578 d1 = tj - old_knots->GetKnot(i);
00579 d2 = old_knots->GetKnot(i + p + ord - v - 1) - tj;
00580
00581 beta = ah->GetKnot(AhIndex(p - 1, i + ord - muprim - 1,ord)) /
00582 (d1 + d2);
00583
00584 ah->PutKnot(AhIndex(p, i + ord - muprim - 2,ord), d2 * beta + beta1) ;
00585 beta1 = d1 * beta;
00586 }
00587
00588 ah->PutKnot(AhIndex(p, iu + ord - muprim - 1,ord), beta1);
00589
00590 if (iu < muprim)
00591 {
00592 register G4double kkk;
00593 register G4double ahv;
00594
00595 kkk = old_knots->GetKnot(n1 - 1 + ord);
00596 ahv = AhIndex (p - 1, iu + ord - muprim,ord);
00597 ah->PutKnot(AhIndex(p, iu + ord - muprim - 1,ord),
00598 beta1 + (kkk - tj) * ahv /
00599 (kkk - old_knots->GetKnot(iu + 1)));
00600 }
00601 }
00602
00603 delete oslo_m->GetKnotVector();
00604 oslo_m->SetKnotVector(new G4KnotVector(v+1));
00605 oslo_m->SetOffset(Amax(muprim - v, 0));
00606 oslo_m->SetSize(v);
00607
00608 for ( i = v, p = 0; i >= 0; i--)
00609 oslo_m->GetKnotVector()
00610 ->PutKnot( p++, ah->GetKnot(AhIndex (v, (ord-1) - i,ord)) );
00611
00612 }
00613
00614 delete ah;
00615 delete newknots;
00616 oslo_m->SetNextNode(oslo_m);
00617 oslo_m = o_ptr;
00618 }
00619
00620 void G4ProjectedSurface::MapSurface(G4ProjectedSurface* srf)
00621 {
00622
00623
00624
00625
00626 register G4ControlPoints *c_ptr;
00627 register G4OsloMatrix *o_ptr;
00628 register G4ControlPoints* new_pts;
00629 register G4ControlPoints* old_pts;
00630
00631 new_pts = srf->ctl_points;
00632
00633
00634
00635
00636
00637 old_pts = new G4ControlPoints(*ctl_points);
00638 register G4int j,
00639 i;
00640 c_ptr = new_pts;
00641
00642 register G4int size;
00643
00644
00645 if(!dir)
00646 size=new_pts->GetRows();
00647 else
00648 size=new_pts->GetCols();
00649
00650 for( register G4int a=0; a<size;a++)
00651 {
00652 if ( lower != 0)
00653 for ( i = 0, o_ptr = oslo_m; i < lower; i++, o_ptr = o_ptr->GetNextNode()){;}
00654 else
00655 o_ptr = oslo_m;
00656
00657 if(!dir)
00658 {
00659 for ( j = lower; j < upper; j++, o_ptr = o_ptr->GetNextNode())
00660 {
00661 register G4double o_scale;
00662 register G4int x;
00663 x=a;
00664
00665
00666
00667
00668
00669 G4Point3D o_pts = old_pts->Get3D(x, o_ptr->GetOffset());
00670 G4Point3D tempc = c_ptr->Get3D(j/upper, (j)%upper-lower);
00671 o_scale = o_ptr->GetKnotVector()->GetKnot(0);
00672
00673 tempc.setX(o_pts.x() * o_scale);
00674 tempc.setY(o_pts.y() * o_scale);
00675
00676 for ( i = 1; i <= o_ptr->GetSize(); i++)
00677 {
00678 o_scale = o_ptr->GetKnotVector()->GetKnot(i);
00679
00680
00681
00682
00683
00684
00685
00686 o_pts = old_pts->Get3D(x,i+o_ptr->GetOffset());
00687 tempc.setX(tempc.x() + o_scale * o_pts.x());
00688 tempc.setY(tempc.y() + o_scale * o_pts.y());
00689 }
00690
00691 c_ptr->put(a,(j)%upper-lower,tempc);
00692 }
00693 }
00694 else
00695 {
00696 for ( j = lower; j < upper; j++, o_ptr = o_ptr->GetNextNode())
00697 {
00698 register G4double o_scale;
00699 register G4int x;
00700 x=a;
00701
00702
00703
00704
00705
00706 G4Point3D o_pts = old_pts->Get3D(o_ptr->GetOffset(),x);
00707 G4Point3D tempc = c_ptr->Get3D((j)%upper-lower, j/upper);
00708
00709 o_scale = o_ptr->GetKnotVector()->GetKnot(0);
00710
00711 tempc.setX(o_pts.x() * o_scale);
00712 tempc.setY(o_pts.y() * o_scale);
00713
00714 for ( i = 1; i <= o_ptr->GetSize(); i++)
00715 {
00716 o_scale = o_ptr->GetKnotVector()->GetKnot(i);
00717 o_pts= old_pts->Get3D(i+o_ptr->GetOffset(),a);
00718
00719 tempc.setX(tempc.x() + o_scale * o_pts.x());
00720 tempc.setY(tempc.y() + o_scale * o_pts.y());
00721 }
00722
00723 c_ptr->put((j)%upper-lower,a,tempc);
00724 }
00725 }
00726 }
00727
00728 delete old_pts;
00729 }