G4ProjectedSurface.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // ----------------------------------------------------------------------
00030 // GEANT 4 class source file
00031 //
00032 // G4ProjectedSurface.cc
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   // Copies the projected surface into a bezier surface
00084   // and adds it to the List of bezier surfaces.
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   // Finds the bounds of the 2D-projected nurb iow
00103   // calculates the bounds for a bounding rectangle
00104   // to the surface. The bounding rectangle is used
00105   // for a preliminary check of intersection.
00106   
00107   // Loop to search the whole control point mesh
00108   // for the minimum and maximum values for x and y.
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 /* L. Broglia
00123       G4Point2d& tmp = (G4Point2d&)ctl_points->get(a,b);
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   // Check wether the surface is a bezier surface by checking
00149   // if internal knots exist.
00150   if(CheckBezier())
00151   {
00152     // Make it a G4BezierSurface -object and add it to the bezier
00153     // surface List     
00154     CopySurface();  
00155         
00156     // Retrieve a pointer to the newly added surface iow the
00157     // last in the List
00158     G4BezierSurface* bez_ptr = (G4BezierSurface*)bezier_list->GetLastSurface();
00159     
00160     // Do the first clip to the bezier.
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         // The clipping routine selected a larger Area than one
00171         // knot interval which indicates that we have a case of
00172         // multiple intersections. The projected surface has to
00173         // be split again in order to separate the intersections
00174         // to different surfaces.
00175 
00176         // Check tolerance of clipping
00177         //          G4cout << "\nClip Area too big -> Split";
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           // The ray intersects with the bounding box
00189           // but not with the surface itself.  
00190           //            G4cout << "\nConvex hull missed.";
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     // Other G4ThreeVec clipping and testing.
00206     bez_ptr->ChangeDir();//bez->dir = !bez_ptr->dir;
00207     bez_ptr->ClipSurface();
00208     //  G4cout<<"\nSMIN: " <<  bez_ptr->smin << "  SMAX: " 
00209     //        <<  bez_ptr->smax << " DIR: " << bez_ptr->dir;
00210            
00211     dMin = bez_ptr->SMin();
00212     dMax = bez_ptr->SMax();
00213     dMaxMinusdMin = dMax-dMin;  
00214     
00215     if((dMaxMinusdMin > kCarTolerance ))// ||
00216       //           (dMaxMinusdMin < -kCarTolerance))
00217     {
00218       if( (dMaxMinusdMin) > 0.8 )
00219       {
00220         //          G4cout << "\nClip Area too big -> Split";       
00221         dir = bez_ptr->dir;//1.2 klo 18.30
00222         
00223         //          dir=!dir;
00224         bezier_list->RemoveSurface(bez_ptr);
00225         SplitNURBSurface();
00226         return;
00227                  //}
00228       }
00229       else
00230         if( dMin > 1.0 || dMax < 0.0 )
00231         {
00232           //            G4cout << "\nConvex hull missed.";
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     // Split the surface into two new surfaces. The G4ThreeVec
00251     // is set in the CheckBezier function. 
00252     //  G4cout << "\nNot a bezier surface -> Split";            
00253     SplitNURBSurface();
00254   }
00255 }
00256 
00257 
00258 G4int G4ProjectedSurface::CheckBezier()
00259 {
00260   // Checks if the surface is a bezier surface by
00261   // checking wether internal knots exist. If no internal
00262   // knots exist the quantity of knots is 2*order of the
00263   // surface. Returns 1 if the surface 
00264   // is a bezier.
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   // Divides the surface in two parts. Uses the oslo-algorithm to calculate
00279   // the new knotvectors and controlpoints for  the subsurfaces.
00280 
00281   //    G4cout << "\nProjected splitted.";
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     //srf1->dir=ROW;
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     //srf2->dir = ROW;
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 // G4ThreeVec = col
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     //srf1->dir = COL;
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     //srf2->dir = COL;
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   // Check that surfaces are ok.
00406   G4int col_size = srf1->ctl_points->GetCols();
00407   G4int row_size = srf1->ctl_points->GetRows();
00408 
00409 /* L. Broglia  
00410   // get three cornerpoints of the controlpoint mesh.
00411   G4Point2d pt1 = srf1->ctl_points->get(0,0);    
00412   G4Point2d pt2 =  srf1->ctl_points->get(0,col_size-1);    
00413   G4Point2d pt3 =  srf1->ctl_points->get(row_size-1,0);    
00414   
00415   // Calc distance between points
00416   G4double pointDist1 = pt1.Distance(pt2);
00417   G4double pointDist2 = pt1.Distance(pt3);
00418 */
00419 
00420   // get three cornerpoints of the controlpoint mesh.
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   // Calc distance squared between points
00426   G4double pointDist1 = pt1.distance2(pt2);
00427   G4double pointDist2 = pt1.distance2(pt3);
00428 
00429 
00430   // Add surfaces to List of projected surfaces    
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 /* L. Broglia      
00440   // get three cornerpoints of the controlpoint mesh.
00441   pt1 = srf2->ctl_points->get(0,0);    
00442   pt2 =  srf2->ctl_points->get(0,col_size-1);    
00443   pt3 =  srf2->ctl_points->get(row_size-1,0);    
00444   
00445   // Calc distance between points
00446   pointDist1 = pt1.Distance(pt2);
00447   pointDist2 = pt1.Distance(pt3);
00448 */
00449 
00450   // get three cornerpoints of the controlpoint mesh.
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   // Calc distance squared between points
00456   pointDist1 = pt1.distance2(pt2);
00457   pointDist2 = pt1.distance2(pt3);
00458 
00459   // Add surfaces to List of projected surfaces    
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   // This algorithm is described in the paper "Making the Oslo-algorithm
00473   // more efficient" in SIAM J.NUMER.ANAL. Vol.23, No. 3, June '86
00474   // Calculates the oslo-matrix , which is used in mapping the new
00475   // knotvector- and controlpoint-values.
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     //while (old_knots->GetKnot(mu + 1) <= new_knots->GetKnot(j))
00527     while ( (new_knots->GetKnot(j) - old_knots->GetKnot(mu + 1)) > 
00528             kCarTolerance                                            )
00529       mu = mu + 1;              // find the bounding mu 
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       // if (new_knots->GetKnot(j + p) == old_knots->GetKnot(ih))
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   // This algorithm is described in the paper "Making the Oslo-algorithm
00623   // more efficient" in SIAM J.NUMER.ANAL. Vol.23, No. 3, June '86
00624   // Maps the new controlpoints into the new surface.
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   // Copy the old points so they can be used in calculating the new ones.
00634   // In this version, where the splitted surfaces are given
00635   // as parameters the copying is not necessary.
00636   
00637   old_pts = new G4ControlPoints(*ctl_points);
00638   register G4int j,    //  j loop 
00639                  i;    //  oslo loop 
00640   c_ptr = new_pts;
00641   
00642   register G4int size; // The number of rows or columns, 
00643                        // depending on processing order
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)// Direction ROW
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 /* L. Broglia   
00666         G4Point2d o_pts = (G4Point2d&)old_pts->get(x,o_ptr->GetOffset());
00667         G4Point2d tempc = (G4Point2d&)c_ptr->get(j/upper,(j)%upper-lower);
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 /* L. Broglia     
00681           o_pts = (G4Point2d&)old_pts->get(x,i+o_ptr->GetOffset());
00682           tempc.X(tempc.X() + o_scale * o_pts.X());
00683           tempc.Y(tempc.Y() + o_scale * o_pts.Y());
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  // dir = COL
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 /* L.Broglia
00703         G4Point2d o_pts = (G4Point2d&)old_pts->get(o_ptr->GetOffset(),x);
00704         G4Point2d tempc = (G4Point2d&)c_ptr->get((j)%upper-lower,j/upper);
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 }

Generated on Mon May 27 17:49:26 2013 for Geant4 by  doxygen 1.4.7