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
00037
00038
00039
00040
00041 #include "G4BezierSurface.hh"
00042 #include "G4ConvexHull.hh"
00043
00044 G4double G4BezierSurface::Tolerance=0;
00045 G4int G4BezierSurface::Clips=0;
00046 G4int G4BezierSurface::Splits=0;
00047
00048 G4BezierSurface::G4BezierSurface()
00049 : G4Surface(), bezier_list(0), smin(0.), smax(0.),
00050 average_u(0.), average_v(0.), dir(0), u_knots(0), v_knots(0),
00051 ctl_points(0), new_knots(0), ord(0), oslo_m(0), lower(0), upper(0),
00052 u_min(0.), u_max(0.), v_min(0.), v_max(0.), old_points(0)
00053 {
00054 order[0]=0; order[1]=0;
00055 u[0]=0.; u[1]=0.;
00056 v[0]=0.; v[1]=0.;
00057 }
00058
00059 G4BezierSurface::~G4BezierSurface()
00060 {
00061 delete u_knots;
00062 delete v_knots;
00063 delete new_knots;
00064 delete ctl_points;
00065 delete old_points;
00066
00067 G4OsloMatrix* temp_oslo = oslo_m;
00068
00069 while(oslo_m != (G4OsloMatrix*)0)
00070 {
00071 oslo_m = oslo_m->GetNextNode();
00072 delete temp_oslo;
00073 temp_oslo = oslo_m;
00074 }
00075
00076 delete oslo_m;
00077 delete bbox;
00078 }
00079
00080 G4BezierSurface::G4BezierSurface(const G4BezierSurface&)
00081 : G4Surface(), bezier_list(0), smin(0.), smax(0.),
00082 average_u(0.), average_v(0.), dir(0), u_knots(0), v_knots(0),
00083 ctl_points(0), new_knots(0), ord(0), oslo_m(0), lower(0), upper(0),
00084 u_min(0.), u_max(0.), v_min(0.), v_max(0.), old_points(0)
00085 {
00086 order[0]=0; order[1]=0;
00087 u[0]=0.; u[1]=0.;
00088 v[0]=0.; v[1]=0.;
00089 }
00090
00091 G4Vector3D G4BezierSurface::SurfaceNormal(const G4Point3D&) const
00092 {
00093 return G4Vector3D(0,0,0);
00094 }
00095
00096 G4int G4BezierSurface::ClipBothDirs()
00097 {
00098 dir = ROW;
00099 ClipSurface();
00100
00101
00102
00103 if(smin > 1.0 || smax < 0.0)
00104 {
00105 bezier_list->RemoveSurface(this);
00106 return 1;
00107 }
00108 else
00109 if((smax - smin) > 0.8)
00110 {
00111 SplitNURBSurface();
00112 return 0;
00113 }
00114
00115 LocalizeClipValues();
00116 SetValues();
00117
00118
00119 dir = COL;
00120 ClipSurface();
00121
00122
00123 if(smin > 1.0 || smax < 0.0)
00124 {
00125 bezier_list->RemoveSurface(this);
00126 return 1;
00127 }
00128 else
00129 if((smax - smin) > 0.8)
00130 {
00131 SplitNURBSurface();
00132 return 0;
00133 }
00134
00135 LocalizeClipValues();
00136 SetValues();
00137 CalcAverage();
00138 return 1;
00139 }
00140
00141
00142 void G4BezierSurface::CalcBBox()
00143 {
00144
00145
00146
00147
00148 G4Point3D box_min = G4Point3D(PINFINITY);
00149 G4Point3D box_max = G4Point3D(-PINFINITY);
00150
00151
00152
00153
00154 for(register G4int a = ctl_points->GetRows()-1; a>=0;a--)
00155 for(register G4int b = ctl_points->GetCols()-1; b>=0;b--)
00156 {
00157
00158
00159
00160
00161
00162
00163
00164 G4Point3D tmp = ctl_points->Get3D(a,b);
00165 if((box_min.x()) > (tmp.x())) box_min.setX(tmp.x());
00166 if((box_max.x()) < (tmp.x())) box_max.setX(tmp.x());
00167 if((box_min.y()) > (tmp.y())) box_min.setY(tmp.y());
00168 if((box_max.y()) < (tmp.y())) box_max.setY(tmp.y());
00169 }
00170
00171 bbox = new G4BoundingBox3D(box_min, box_max);
00172 }
00173
00174
00175 void G4BezierSurface::CalcAverage()
00176 {
00177
00178 average_u = (u_min + u_max)/2.0;
00179 average_v = (v_min + v_max)/2.0;
00180 }
00181
00182
00183 void G4BezierSurface::CalcDistance(const G4Point3D& ray_start)
00184 {
00185
00186
00187 distance = ((((ray_start.x() - average_pt.x())*
00188 (ray_start.x() - average_pt.x()))+
00189 ((ray_start.y() - average_pt.y())*
00190 (ray_start.y() - average_pt.y()))+
00191 ((ray_start.z() - average_pt.z())*
00192 (ray_start.z() - average_pt.z()))));
00193 }
00194
00195
00196 void G4BezierSurface::SetValues()
00197 {
00198 if(dir)
00199 {
00200 v_min = smin;
00201 v_max = smax;
00202 }
00203 else
00204 {
00205 u_min = smin;
00206 u_max = smax;
00207 }
00208 }
00209
00210
00211 G4int G4BezierSurface::BIntersect(G4SurfaceList& bez_list)
00212 {
00213 bezier_list = &bez_list;
00214 G4int clip_regions = 0;
00215
00216 do
00217 {
00218
00219 CalcBBox();
00220
00221
00222
00223
00224
00225
00226
00227
00228 if(!bbox->GetTestResult())
00229 return 0;
00230
00231
00232
00233
00234
00235
00236 GetClippedRegionFromSurface();
00237 clip_regions++;
00238
00239
00240
00241 RefineSurface();
00242
00243
00244 u_min = u_knots->GetKnot(0);
00245 u_max = u_knots->GetKnot(u_knots->GetSize() - 1);
00246 v_min = v_knots->GetKnot(0);
00247 v_max = v_knots->GetKnot(v_knots->GetSize() - 1);
00248
00249
00250
00251 if( (u_max - u_min) < (v_max - v_min) )
00252 dir = 1;
00253 else
00254 dir = 0;
00255
00256
00257 ClipSurface();
00258
00259
00260
00261
00262 if( smin > 1.0 || smax < 0.0 )
00263 {
00264
00265
00266 return 0;
00267 }
00268
00269 if( (smax - smin) > 0.8)
00270 {
00271
00272
00273 SplitNURBSurface();
00274
00275
00276
00277
00278
00279
00280
00281
00282 return 2;
00283 }
00284
00285
00286
00287 LocalizeClipValues();
00288
00289
00290
00291 } while ((u_max - u_min > Tolerance) || (v_max - v_min) > Tolerance);
00292
00293 SetValues();
00294
00295
00296
00297 return 1;
00298 }
00299
00300
00301 void G4BezierSurface::ClipSurface()
00302 {
00303
00304
00305
00306
00307
00308
00309
00310 register G4int i,j;
00311 register G4int col_size = ctl_points->GetCols();
00312 register G4int row_size = ctl_points->GetRows();
00313
00314 G4ConvexHull *ch_tmp= new G4ConvexHull(0,1.0e8,-1.0e8);
00315 G4ConvexHull *ch_ptr=0, *ch_first=0;
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326 G4Point3D pt1 = ctl_points->Get3D(0,0);
00327 G4Point3D pt2 = ctl_points->Get3D(0,col_size-1);
00328 G4Point3D pt3 = ctl_points->Get3D(row_size-1,0);
00329 G4Point3D pt4 = ctl_points->Get3D(row_size-1,col_size-1);
00330 G4Point3D v1,v2,v3;
00331
00332 if ( dir == ROW)
00333 {
00334
00335 v1 = (pt1 - pt3);
00336
00337
00338 v2 = (pt2 - pt4);
00339
00340
00341 }
00342 else
00343 {
00344 v1 = pt1 - pt2;
00345 v2 = pt3 - pt4;
00346
00347
00348
00349
00350 }
00351
00352
00353
00354
00355 v3 = v1 + v2 ;
00356
00357 smin = 1.0e8;
00358 smax = -1.0e8;
00359
00360 G4double norm = std::sqrt(v3.x() * v3.x() + v3.y() * v3.y());
00361 if(!norm)
00362 {
00363 G4cout << "\nNormal zero!";
00364 G4cout << "\nLINE & DIR: " << line.x() << " " << line.y() << " " << dir;
00365 G4cout << "\n";
00366
00367 if((std::abs(line.x())) > kCarTolerance)
00368 line.setX(-line.x());
00369 else
00370 if((std::abs(line.y())) > kCarTolerance)
00371 line.setY(-line.y());
00372 else
00373 {
00374 G4cout << "\n RETURNING FROm CLIP..";
00375 smin = 0; smax = 1; delete ch_tmp;
00376 return;
00377 }
00378
00379 G4cout << "\nCHANGED LINE & DIR: " << line.x() << " "
00380 << line.y() << " " << dir;
00381 }
00382 else
00383 {
00384 line.setX( v3.y() / norm);
00385 line.setY(-v3.x() / norm);
00386 }
00387
00388
00389
00390
00391
00392
00393 if( dir == ROW)
00394 {
00395
00396 for(G4int a = 0; a < col_size; a++)
00397 {
00398 ch_ptr = new G4ConvexHull(a/(col_size - 1.0),1.0e8,-1.0e8);
00399 if(! a)
00400 {
00401 ch_first=ch_ptr; delete ch_tmp;
00402 }
00403 else ch_tmp->SetNextHull(ch_ptr);
00404
00405 ch_tmp=ch_ptr;
00406 }
00407
00408 ch_ptr=ch_first;
00409 register G4double value;
00410
00411
00412
00413
00414 for( G4int h = 0; h < row_size; h++)
00415 {
00416 for(G4int k = 0; k < col_size; k++)
00417 {
00418
00419
00420
00421
00422 G4Point3D coordstmp = ctl_points->Get3D(h,k);
00423 value = - ((coordstmp.x() * line.x() + coordstmp.y() * line.y()));
00424
00425 if( value <= (ch_ptr->GetMin()+kCarTolerance)) ch_ptr->SetMin(value);
00426 if( value >= (ch_ptr->GetMax()-kCarTolerance)) ch_ptr->SetMax(value);
00427
00428 ch_ptr=ch_ptr->GetNextHull();
00429 }
00430
00431 ch_ptr=ch_first;
00432 }
00433
00434 ch_ptr=ch_first;
00435
00436
00437
00438
00439
00440 for(G4int l = 0; l < col_size - 1; l++)
00441 {
00442 ch_tmp=ch_ptr->GetNextHull();
00443 for(G4int q = l+1; q < col_size; q++)
00444 {
00445 register G4double d, param1, param2;
00446 param1 = ch_ptr->GetParam();
00447 param2 = ch_tmp->GetParam();
00448
00449 if(ch_tmp->GetMax() - ch_ptr->GetMax())
00450 {
00451 d = Findzero( param1, param2, ch_ptr->GetMax(), ch_tmp->GetMax());
00452 if( d <= (smin + kCarTolerance) ) smin = d * .99;
00453 if( d >= (smax - kCarTolerance) ) smax = d * .99 + .01;
00454 }
00455
00456 if(ch_tmp->GetMin() - ch_ptr->GetMin())
00457 {
00458 d = Findzero( param1, param2, ch_ptr->GetMin(), ch_tmp->GetMin());
00459 if( d <= (smin + kCarTolerance)) smin = d * .99;
00460 if( d >= (smax - kCarTolerance)) smax = d * .99 + .01;
00461 }
00462
00463 ch_tmp=ch_tmp->GetNextHull();
00464 }
00465
00466 ch_ptr=ch_ptr->GetNextHull();
00467 }
00468
00469 ch_ptr=ch_first;
00470
00471 if (smin <= 0.0) smin = 0.0;
00472 if (smax >= 1.0) smax = 1.0;
00473
00474 if ( (ch_ptr)
00475 && (Sign(ch_ptr->GetMin()) != Sign(ch_ptr->GetMax()))) smin = 0.0;
00476
00477 i = Sign(ch_tmp->GetMin());
00478 j = Sign(ch_tmp->GetMax());
00479
00480 if ( std::abs(i-j) > kCarTolerance ) smax = 1.0;
00481
00482
00483 }
00484 else
00485 {
00486 for(G4int n = 0; n < row_size; n++)
00487 {
00488 ch_ptr = new G4ConvexHull(n/(row_size - 1.0),1.0e8,-1.0e8);
00489 if(!n)
00490 {
00491 ch_first=ch_ptr; delete ch_tmp;
00492 }
00493 else ch_tmp->SetNextHull(ch_ptr);
00494
00495 ch_tmp=ch_ptr;
00496 }
00497
00498 ch_ptr=ch_first;
00499
00500 for( G4int o = 0; o < col_size; o++)
00501 {
00502 for(G4int p = 0; p < row_size; p++)
00503 {
00504 register G4double value;
00505
00506
00507
00508
00509
00510 G4Point3D coordstmp = ctl_points->Get3D(p,o);
00511 value = - ((coordstmp.x() * line.x() + coordstmp.y() * line.y()));
00512
00513 if( value <= (ch_ptr->GetMin()+kCarTolerance)) ch_ptr->SetMin(value);
00514 if( value >= (ch_ptr->GetMax()-kCarTolerance)) ch_ptr->SetMax(value);
00515
00516 ch_ptr=ch_ptr->GetNextHull();
00517 }
00518
00519 ch_ptr=ch_first;
00520 }
00521
00522 ch_ptr=ch_first;
00523 delete ch_tmp;
00524
00525 for(G4int q = 0; q < row_size - 1; q++)
00526 {
00527 ch_tmp=ch_ptr->GetNextHull();
00528 for(G4int r = q+1; r < row_size; r++)
00529 {
00530 register G4double param1 = ch_ptr->GetParam();
00531 register G4double param2 = ch_tmp->GetParam();
00532 register G4double d;
00533
00534 if(ch_tmp->GetMax() - ch_ptr->GetMax())
00535 {
00536 d = Findzero( param1, param2, ch_ptr->GetMax(), ch_tmp->GetMax());
00537 if( d <= (smin + kCarTolerance) ) smin = d * .99;
00538 if( d >= (smax - kCarTolerance) ) smax = d * .99 + .01;
00539 }
00540
00541 if(ch_tmp->GetMin()-ch_ptr->GetMin())
00542 {
00543 d = Findzero( param1, param2, ch_ptr->GetMin(), ch_tmp->GetMin());
00544 if( d <= (smin + kCarTolerance) ) smin = d * .99;
00545 if( d >= (smax - kCarTolerance) ) smax = d * .99 + .01;
00546 }
00547
00548 ch_tmp=ch_tmp->GetNextHull();
00549 }
00550
00551 ch_ptr=ch_ptr->GetNextHull();
00552 }
00553
00554 ch_tmp=ch_ptr;
00555 ch_ptr=ch_first;
00556
00557 if (smin <= 0.0) smin = 0.0;
00558 if (smax >= 1.0) smax = 1.0;
00559
00560 if ( (ch_ptr)
00561 && (Sign(ch_ptr->GetMin()) != Sign(ch_ptr->GetMax()))) smin = 0.0;
00562
00563 if ( ch_tmp )
00564 {
00565 i = Sign(ch_tmp->GetMin());
00566 j = Sign(ch_tmp->GetMax());
00567 if ( (std::abs(i-j) > kCarTolerance)) smax = 1.0;
00568 }
00569 }
00570
00571 ch_ptr=ch_first;
00572 while(ch_ptr && (ch_ptr!=ch_ptr->GetNextHull()))
00573 {
00574 ch_tmp=ch_ptr;
00575 ch_ptr=ch_ptr->GetNextHull();
00576 delete ch_tmp;
00577 }
00578
00579 delete ch_ptr;
00580
00581
00582 Clips++;
00583 }
00584
00585
00586 void G4BezierSurface::GetClippedRegionFromSurface()
00587 {
00588
00589
00590
00591
00592
00593
00594 delete new_knots;
00595 if ( dir == ROW)
00596 {
00597 new_knots = new G4KnotVector(GetOrder(0) * 2);
00598 for (register G4int i = 0; i < GetOrder(0); i++)
00599 {
00600 new_knots->PutKnot(i, smin);
00601 new_knots->PutKnot(i+ GetOrder(0), smax);
00602 }
00603 }
00604 else
00605 {
00606 new_knots = new G4KnotVector( GetOrder(1) * 2);
00607 for ( register G4int i = 0; i < GetOrder(1); i++)
00608 {
00609 new_knots->PutKnot(i, smin);
00610 new_knots->PutKnot(i+ GetOrder(1), smax);
00611 }
00612 }
00613 }
00614
00615
00616 void G4BezierSurface::RefineSurface()
00617 {
00618
00619
00620
00621 delete old_points;
00622 if (dir == ROW)
00623 {
00624
00625 ord = GetOrder(0);
00626 CalcOsloMatrix();
00627 for(register G4int a=0;a<new_knots->GetSize();a++)
00628 u_knots->PutKnot(a, new_knots->GetKnot(a));
00629
00630 lower = 0;
00631 upper = new_knots->GetSize() - GetOrder(0);
00632
00633
00634 old_points = new G4ControlPoints(*ctl_points);
00635 MapSurface(this);
00636 }
00637 else
00638 {
00639 ord = GetOrder(1);
00640 CalcOsloMatrix ();
00641 for(register G4int a=0;a < new_knots->GetSize();a++)
00642 v_knots->PutKnot(a, new_knots->GetKnot(a));
00643
00644
00645 old_points = new G4ControlPoints(*ctl_points);
00646
00647
00648 register G4int cols = ctl_points->GetCols();
00649 delete ctl_points;
00650
00651 ctl_points = new G4ControlPoints(2,(new_knots->GetSize()-
00652 GetOrder(1)),cols);
00653 lower = 0;
00654 upper = new_knots->GetSize() - GetOrder(1);
00655 MapSurface(this);
00656 }
00657 }
00658
00659
00660 void G4BezierSurface::CalcOsloMatrix()
00661 {
00662
00663
00664
00665
00666
00667 register G4KnotVector *ah;
00668 register G4KnotVector *newknots;
00669 register G4int i;
00670 register G4int j;
00671 register G4int mu, muprim;
00672 register G4int vv, p;
00673 register G4int iu, il, ih, n1;
00674 register G4int ahi;
00675 register G4double beta1;
00676 register G4double tj;
00677
00678 ah = new G4KnotVector(ord*(ord + 1)/2);
00679 newknots = new G4KnotVector(ord * 2 );
00680
00681 n1 = new_knots->GetSize() - ord;
00682 mu = 0;
00683
00684 if(oslo_m!=(G4OsloMatrix*)0)
00685 {
00686 G4OsloMatrix* tmp;
00687
00688
00689 while(oslo_m!=(G4OsloMatrix*)0)
00690 {
00691 tmp=oslo_m->GetNextNode();delete oslo_m; oslo_m=tmp;
00692 }
00693 }
00694
00695 delete oslo_m;
00696 oslo_m = new G4OsloMatrix();
00697
00698 register G4OsloMatrix* o_ptr = oslo_m;
00699
00700 register G4KnotVector* old_knots;
00701 if(dir)
00702 old_knots = v_knots;
00703 else
00704 old_knots = u_knots;
00705
00706 for (j = 0; j < n1; j++)
00707 {
00708 if ( j != 0 )
00709 {
00710 oslo_m->SetNextNode(new G4OsloMatrix());
00711 oslo_m = oslo_m->GetNextNode();
00712 }
00713
00714 while (old_knots->GetKnot(mu + 1) <= new_knots->GetKnot(j))
00715 mu = mu + 1;
00716
00717 i = j + 1;
00718 muprim = mu;
00719
00720 while ((new_knots->GetKnot(i) == old_knots->GetKnot(muprim)) &&
00721 i < (j + ord))
00722 {
00723 i++;
00724 muprim--;
00725 }
00726
00727 ih = muprim + 1;
00728
00729 for (vv = 0, p = 1; p < ord; p++)
00730 {
00731 if (new_knots->GetKnot(j + p) == old_knots->GetKnot(ih))
00732 ih++;
00733 else
00734 newknots->PutKnot(++vv - 1,new_knots->GetKnot(j + p));
00735 }
00736
00737 ahi = AhIndex(0, ord - 1,ord);
00738 ah->PutKnot(ahi, 1.0);
00739
00740 for (p = 1; p <= vv; p++)
00741 {
00742 beta1 = 0.0;
00743 tj = newknots->GetKnot(p-1);
00744
00745 if (p - 1 >= muprim)
00746 {
00747 beta1 = AhIndex(p - 1, ord - muprim,ord);
00748 beta1 = ((tj - old_knots->GetKnot(0)) * beta1) /
00749 (old_knots->GetKnot(p + ord - vv) - old_knots->GetKnot(0));
00750 }
00751
00752 i = muprim - p + 1;
00753 il = Amax (1, i);
00754 i = n1 - 1 + vv - p;
00755 iu = Amin (muprim, i);
00756
00757 for (i = il; i <= iu; i++)
00758 {
00759 register G4double d1, d2;
00760 register G4double beta;
00761
00762 d1 = tj - old_knots->GetKnot(i);
00763 d2 = old_knots->GetKnot(i + p + ord - vv - 1) - tj;
00764
00765 beta = ah->GetKnot(AhIndex(p - 1, i + ord - muprim - 1,ord)) /
00766 (d1 + d2);
00767
00768
00769 ah->PutKnot(AhIndex(p, i + ord - muprim - 2,ord), d2 * beta + beta1) ;
00770 beta1 = d1 * beta;
00771 }
00772
00773 ah->PutKnot(AhIndex(p, iu + ord - muprim - 1,ord), beta1);
00774
00775 if (iu < muprim)
00776 {
00777 register G4double kkk;
00778 register G4double ahv;
00779
00780 kkk = old_knots->GetKnot(n1 - 1 + ord);
00781 ahv = AhIndex (p - 1, iu + ord - muprim,ord);
00782 ah->PutKnot(AhIndex(p, iu + ord - muprim - 1,ord),
00783 beta1 + (kkk - tj) * ahv /
00784 (kkk - old_knots->GetKnot(iu + 1)));
00785 }
00786 }
00787
00788
00789 G4OsloMatrix* temp_oslo = oslo_m;
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804 while(oslo_m != (G4OsloMatrix*)0)
00805 {
00806 oslo_m = oslo_m->GetNextNode();
00807 delete temp_oslo;
00808 temp_oslo = oslo_m;
00809 }
00810
00811 delete oslo_m;
00812
00813
00814 oslo_m = new G4OsloMatrix(vv+1, Amax(muprim - vv,0), vv);
00815
00816 for ( i = vv, p = 0; i >= 0; i--)
00817 oslo_m->GetKnotVector()
00818 ->PutKnot ( p++, ah->GetKnot(AhIndex (vv, (ord-1) - i,ord)));
00819
00820 }
00821
00822 delete ah;
00823 delete newknots;
00824 oslo_m->SetNextNode(0);
00825 oslo_m = o_ptr;
00826 }
00827
00828
00829 void G4BezierSurface::MapSurface(G4Surface*)
00830 {
00831
00832
00833
00834
00835 register G4ControlPoints *c_ptr;
00836 register G4OsloMatrix *o_ptr;
00837 register G4ControlPoints* new_pts;
00838 register G4ControlPoints* old_pts;
00839
00840 new_pts = ctl_points;
00841
00842
00843
00844 old_pts = old_points;
00845 register G4int j,
00846 i;
00847
00848 c_ptr = new_pts;
00849 register G4int size;
00850
00851
00852 if(!dir)
00853 size=new_pts->GetRows();
00854 else
00855 size=new_pts->GetCols();
00856
00857 for(G4int a=0; a<size;a++)
00858 {
00859 if ( lower != 0)
00860 {
00861 for ( i = 0, o_ptr = oslo_m;
00862 i < lower;
00863 i++, o_ptr = o_ptr->GetNextNode()){;}
00864 }
00865 else
00866 {
00867 o_ptr = oslo_m;
00868 }
00869
00870 if(!dir)
00871 {
00872 for ( j = lower; j < upper; j++, o_ptr = o_ptr->GetNextNode())
00873 {
00874 register G4double o_scale;
00875 register G4int x;
00876 x=a;
00877
00878
00879
00880
00881
00882 G4Point3D o_pts = old_pts->Get3D(x, o_ptr->GetOffset());
00883 G4Point3D tempc = c_ptr->Get3D(j/upper, (j)%upper-lower);
00884
00885 o_scale = o_ptr->GetKnotVector()->GetKnot(0);
00886
00887 tempc.setX(o_pts.x() * o_scale);
00888 tempc.setY(o_pts.x() * o_scale);
00889
00890 for ( i = 1; i <= o_ptr->GetSize(); i++)
00891 {
00892 o_scale = o_ptr->GetKnotVector()->GetKnot(i);
00893
00894
00895
00896
00897
00898
00899 o_pts = old_pts->Get3D(x, i+o_ptr->GetOffset());
00900 tempc.setX(tempc.x() + o_scale * o_pts.x());
00901 tempc.setY(tempc.y() + o_scale * o_pts.y());
00902
00903 }
00904
00905 c_ptr->put(a,(j)%upper-lower,tempc);
00906 }
00907 }
00908 else
00909 {
00910 for ( j = lower; j < upper; j++, o_ptr = o_ptr->GetNextNode())
00911 {
00912 register G4double o_scale;
00913 register G4int x;
00914 x=a;
00915
00916
00917
00918
00919
00920 G4Point3D o_pts = old_pts->Get3D(o_ptr->GetOffset(), x);
00921 G4Point3D tempc = c_ptr->Get3D((j)%upper-lower,j/upper);
00922
00923 o_scale = o_ptr->GetKnotVector()->GetKnot(0);
00924
00925 tempc.setX(o_pts.x() * o_scale);
00926 tempc.setY(o_pts.y() * o_scale);
00927
00928 for ( i = 1; i <= o_ptr->GetSize(); i++)
00929 {
00930 o_scale = o_ptr->GetKnotVector()->GetKnot(i);
00931
00932
00933
00934 o_pts= old_pts->Get3D(i+o_ptr->GetOffset(),a);
00935 tempc.setX(tempc.x() + o_scale * o_pts.x());
00936 tempc.setY(tempc.y() + o_scale * o_pts.y());
00937 }
00938
00939 c_ptr->put((j)%upper-lower,a,tempc);
00940 }
00941 }
00942 }
00943 }
00944
00945
00946 void G4BezierSurface::SplitNURBSurface()
00947 {
00948
00949
00950
00951
00952
00953 register G4double value;
00954 register G4int i;
00955 register G4int k_index=0;
00956 G4BezierSurface *srf1, *srf2;
00957 G4int nr,nc;
00958
00959 if ( dir == ROW )
00960 {
00961 value = u_knots->GetKnot((u_knots->GetSize()-1)/2);
00962
00963 for( i = 0; i < u_knots->GetSize(); i++)
00964 if( value == u_knots->GetKnot(i) )
00965 {
00966 k_index = i;
00967 break;
00968 }
00969
00970 if ( k_index == 0)
00971 {
00972 value = ( value + u_knots->GetKnot(u_knots->GetSize() -1))/2.0;
00973 k_index = GetOrder(ROW);
00974 }
00975
00976 new_knots = u_knots->MultiplyKnotVector(GetOrder(ROW), value);
00977
00978 ord = GetOrder(ROW);
00979 CalcOsloMatrix();
00980
00981 srf1 = new G4BezierSurface(*this);
00982
00983 srf1->dir=COL;
00984
00985 new_knots->ExtractKnotVector(srf1->u_knots, k_index +
00986 srf1->GetOrder(ROW),0);
00987
00988 nr= srf1->v_knots->GetSize() - srf1->GetOrder(COL);
00989 nc= srf1->u_knots->GetSize() - srf1->GetOrder(ROW);
00990 delete srf1->ctl_points;
00991
00992 srf1->ctl_points= new G4ControlPoints(2, nr, nc);
00993 srf2 = new G4BezierSurface(*this);
00994
00995
00996 srf2->dir = COL;
00997
00998 new_knots->ExtractKnotVector(srf2->u_knots,
00999 new_knots->GetSize(), k_index);
01000
01001 nr= srf2->v_knots->GetSize() - srf2->GetOrder(COL);
01002 nc= srf2->u_knots->GetSize() - srf2->GetOrder(ROW);
01003
01004 delete srf2->ctl_points;
01005 srf2->ctl_points = new G4ControlPoints(2, nr, nc);
01006
01007 lower = 0;
01008 upper = k_index;
01009 MapSurface(srf1);
01010
01011 lower = k_index;
01012 upper = new_knots->GetSize() - srf2->GetOrder(ROW);
01013 MapSurface(srf2);
01014 }
01015 else
01016 {
01017 value = v_knots->GetKnot((v_knots->GetSize() -1)/2);
01018
01019 for( i = 0; i < v_knots->GetSize(); i++)
01020 if( value == v_knots->GetKnot(i))
01021 {
01022 k_index = i;
01023 break;
01024 }
01025 if ( k_index == 0)
01026 {
01027 value = ( value + v_knots->GetKnot(v_knots->GetSize() -1))/2.0;
01028 k_index = GetOrder(COL);
01029 }
01030
01031 new_knots = v_knots->MultiplyKnotVector( GetOrder(COL), value );
01032 ord = GetOrder(COL);
01033
01034 CalcOsloMatrix();
01035
01036 srf1 = new G4BezierSurface(*this);
01037
01038 srf1->dir = ROW;
01039
01040 new_knots->ExtractKnotVector(srf1->v_knots,
01041 k_index + srf1->GetOrder(COL), 0);
01042
01043 nr = srf1->v_knots->GetSize() - srf1->GetOrder(COL);
01044 nc = srf1->u_knots->GetSize() - srf1->GetOrder(ROW);
01045
01046 delete srf1->ctl_points;
01047 srf1->ctl_points = new G4ControlPoints(2, nr, nc);
01048
01049 srf2 = new G4BezierSurface(*this);
01050
01051 srf2->dir = ROW;
01052
01053 new_knots->ExtractKnotVector(srf2->v_knots, new_knots->GetSize(), k_index);
01054
01055 nr = srf2->v_knots->GetSize() - srf2->GetOrder(COL);
01056 nc = srf2->u_knots->GetSize() - srf2->GetOrder(ROW);
01057
01058 delete srf2->ctl_points;
01059 srf2->ctl_points = new G4ControlPoints(2,nr, nc);
01060
01061 lower = 0;
01062 upper = k_index;
01063 MapSurface(srf1);
01064
01065
01066 lower = k_index;
01067 upper = new_knots->GetSize() - srf2->GetOrder(COL);
01068 MapSurface(srf2);
01069 }
01070
01071 bezier_list->AddSurface(srf1);
01072 bezier_list->AddSurface(srf2);
01073 delete new_knots;
01074
01075
01076 Splits++;
01077 }