BFLVorOperator.cxx

Go to the documentation of this file.

//                                                                   //
//authors: Costas Andreopoulos, Elias Athanasopoulos, George Tzanakos//
//                Physics Department, University of Athens           //
//                                                  October 28, 1999 //

#include <fstream>
#include "TObjArray.h"
#include "TRandom.h"
#include "TMatrix.h"
#include "TMath.h"
#include "TVector2.h"

#include "BField/TIntList.h"
//#include "BField/BFLNode.h"
//#include "BField/BFLVtx.h"
//#include "BField/BFLEdge.h"
//#include "BField/BFLPolyg.h"
#include "BField/BFLWingedEdge.h"
#include "BField/BFLVorOperator.h"
#include "BField/BFLVoronoiMaker.h"
   
ClassImp(BFLVorOperator)


// CONSTRUCTOR & DESTRUCTOR
//

BFLVorOperator::BFLVorOperator(BFLWingedEdge * voronoi)
{
// Allocate some memory for this instance variables
  fVor = (BFLWingedEdge *) voronoi;
  fInVtxCache = new TIntList();
  fOutVtxCache = new TIntList(); 
}

BFLVorOperator::~BFLVorOperator()
{
  delete fInVtxCache;
  delete fOutVtxCache;
}

TIntList * BFLVorOperator::RetrieveEdgesIncidentToVtx(Int_t vtx) const
{
// Input  : vertex id
// Output : a pointer to a TList that contains the edges incident
//          to the vertex ( in counterclockwise order )
//
  Int_t k, kstart;
  TIntList * IncidentEdges = new TIntList();

  k = fVor->GetEdgeAroundVtx(vtx);
  kstart = k;

  do {

    IncidentEdges->AddLast(k);
    if( vtx == fVor->GetStartVtx(k) ) {
           k = fVor->GetCcwPredecessor(k);
        } else {
           k = fVor->GetCcwSuccessor(k);
        }

  } while ( k != kstart );

  return IncidentEdges;
}

TIntList * BFLVorOperator::RetrievePolygonsIncidentToVtx(Int_t vtx) const
{
// Input  : vertex id
// Output : a pointer to a TList that contains the polygons
//          incident to the vertex ( in counterclockwise order )
//
  Int_t k, kstart, poly;
  TIntList * IncidentPolygons = new TIntList();

  k = fVor->GetEdgeAroundVtx(vtx);  
  kstart = k;

  do {

    if( vtx == fVor->GetStartVtx(k) ) {
           poly = fVor->GetLeftPolyg(k);
           k = fVor->GetCcwPredecessor(k);
    } else {
           poly = fVor->GetRightPolyg(k);
           k = fVor->GetCcwSuccessor(k);
    }
    IncidentPolygons -> AddLast(poly);

  } while( k != kstart);

  return IncidentPolygons;
}

TIntList * BFLVorOperator::RetrieveEdgesSurrPolygon(Int_t poly) const
{
  Int_t k, kstart;
  TIntList * SurroundingEdges = new TIntList();

  k = fVor->GetEdgeAroundPolyg(poly);
  kstart = k;

  do {

    SurroundingEdges->AddLast(k);

    if(poly == fVor->GetLeftPolyg(k)){
       k = fVor->GetCwSuccessor(k);
     } else {
       k = fVor->GetCwPredecessor(k);
     } 
  } while( k!=kstart );

  return SurroundingEdges;
}


TIntList * BFLVorOperator::RetrieveVtxSurrPolygon(Int_t poly) const
{
  Int_t k, kstart;
  TIntList * SurrVertices = new TIntList();

  k = fVor->GetEdgeAroundPolyg(poly);
  kstart = k;

  do {

    if(poly == fVor->GetLeftPolyg(k)){
       SurrVertices->AddLast(fVor->GetEndVtx(k));
       k = fVor->GetCwSuccessor(k);
    } else {
       SurrVertices->AddLast(fVor->GetStartVtx(k));
       k = fVor->GetCwPredecessor(k);   
    }

  } while( k!=kstart );

  return SurrVertices;
}


Int_t BFLVorOperator::GetFirstIntersectEdge(void) const
{
// Returns the first edge of the polygon: fCurrentPolygID that intersects
// the new Voronoi edge in this polygon. It initiates the boundary 
// growing procedure and is also used to terminate it ( when we return
// back to the edge we had begun from ).
//
  Int_t FirstIntersectEdgeID = 0;
  TIntList * IntersectEdges = new TIntList();

  // Get a list of current polygon's edges 
  TIntList * PolyEdgeIDs = RetrieveEdgesSurrPolygon(fCurrentPolygID);  

  // Select the (2) edges that intersect with the new Voronoi edge 
  // at this region. 
  for(Int_t i=0; i<PolyEdgeIDs->NumberOfElements(); i++){
    Int_t iedge = PolyEdgeIDs->At(i); 
    if(EdgeHasNewVtx(iedge)) IntersectEdges->Add(iedge);
  }
  // check possible fealure (there MUST be 2 intersections)
  if(IntersectEdges->NumberOfElements() != 2) {
     delete PolyEdgeIDs;
     delete IntersectEdges;
     return -1;
  }
 
  // Select the intersect edge that will allow the boundary growing
  // procedure to be performed counterclockwisely with respect to  
  // the newly added generator.

  BFLNode VtxOnFirstEdge = FindNewVtx(IntersectEdges->At(0));
  BFLNode VtxOnSecondEdge = FindNewVtx(IntersectEdges->At(1));  

  if(Clockwise(fCurrentGen,&VtxOnFirstEdge,&VtxOnSecondEdge) == -1) {
     FirstIntersectEdgeID = IntersectEdges->At(0);
  } else FirstIntersectEdgeID = IntersectEdges->At(1);
                                             
  // Explicitly free some allocated memory.  
  delete PolyEdgeIDs;
  delete IntersectEdges;
  
  return FirstIntersectEdgeID;
}

Int_t BFLVorOperator::Clockwise(BFLNode * nA, BFLNode * nB, 
                                BFLNode * nC) const
{
// Returns 1 if going from node nA -> nB -> nC is clockwise
// and -1 if its counterclockwise.
// Special cases : 
//  1 if nB between nA,nC 
// -1 if nA between nB,nC
//  0 if nC between nA,nB
// --------- Adapted from "Algorithms", R.Hedgewick

  Int_t cwise;
  Float_t dxa, dya, dxb, dyb; 

  dxa = ( nB->GetX() ) - ( nA->GetX() );
  dxb = ( nC->GetX() ) - ( nA->GetX() );
  dya = ( nB->GetY() ) - ( nA->GetY() );
  dyb = ( nC->GetY() ) - ( nA->GetY() );

  if (dxa*dyb > dya*dxb) cwise = -1; 
  else if (dxa*dyb < dya*dxb) cwise = 1; 
  else { /* dxa*dyb == dya*dxb */
    if (dxa*dxb < 0 || dya*dyb < 0 ) cwise = -1;
        else if (dxa*dxa + dya*dya >= dxb*dxb + dyb*dyb) cwise = 0; 
        else cwise = 1;
  } 
  
  return cwise;
}

Int_t BFLVorOperator::GetNextIntersectEdge(Int_t PrEdgeID) const
{
// Returns the edge ( other than edge: PrEdgeID) that intersects with 
// the new Voronoi edge on the polygon: fCurrentPolygID. The boundary
// growing procedure is based on successive calls of this methods
// until we return to the edge we had begun from ( and hence the new
// - closed - Voronoi region is formed ).
//
  Int_t NextEdgeID = 0, iedge;
  TIntList * PolyEdgeIDs;

  // Get a list of current polygon's edges
  PolyEdgeIDs = RetrieveEdgesSurrPolygon(fCurrentPolygID);  

  // Remove the previous intersect edge from the list   
  PolyEdgeIDs->Remove(PrEdgeID);

  // Get the intersect edge ( other than PrEdgeID ).
  for(Int_t i = 0; i < PolyEdgeIDs->NumberOfElements(); i++){
    iedge = PolyEdgeIDs->At(i);
    if(EdgeHasNewVtx(iedge)) NextEdgeID = iedge;
  }

  delete PolyEdgeIDs;

  return NextEdgeID;
}

Bool_t BFLVorOperator::EdgeHasNewVtx(Int_t EdgeID) const
{
// Returns TRUE if the edge: EdgeID has a vtx of the region to be
// assigned to the new generator: igen (and hence has to be resized).
// It happens when one of the vertices in within the new region 
// while the other is outside.
//
  Int_t SVtxID, EVtxID;
  Bool_t SVtxIsInside, EVtxIsInside, HasNewVtx;

  // Get start/end vertex IDs of the edge.
  SVtxID = fVor->GetStartVtx(EdgeID);
  EVtxID = fVor->GetEndVtx(EdgeID);

  // Check whether the above vertices are within the region
  // to be assigned to generator : igen. 
  SVtxIsInside = VtxIsInsideNewPolyg(SVtxID); 
  EVtxIsInside = VtxIsInsideNewPolyg(EVtxID);  
   
  // Decide whether edge: EdgeID is to be resized.
  if((SVtxIsInside && !EVtxIsInside)||(!SVtxIsInside && EVtxIsInside)){
     HasNewVtx = kTRUE;
  } else {
     HasNewVtx = kFALSE;
  }
  
  return HasNewVtx;
}

Bool_t BFLVorOperator::EdgeIsInsideNewPolygon(Int_t EdgeID) const
{
// Returns TRUE if the edge: EdgeID has both of its vertices within
// the region to be assigned to the new generator: igen ( and hence
// belongs to the substructure to be deleted ).
//
  Int_t SVtxID, EVtxID;
  Bool_t SVtxIsInside, EVtxIsInside, IsInside;
  
  // Get start/end vertex IDs of the edge.
  SVtxID = fVor->GetStartVtx(EdgeID);
  EVtxID = fVor->GetEndVtx(EdgeID);
  
  // Check whether the above vertices are within the region
  // to be assigned to generator : igen.
  SVtxIsInside = VtxIsInsideNewPolyg(SVtxID);
  EVtxIsInside = VtxIsInsideNewPolyg(EVtxID);  
    
  // Decide whether edge: EdgeID belongs to the substructure
  // to be deleted.
  if(SVtxIsInside && EVtxIsInside) IsInside = kTRUE;
  else IsInside = kFALSE;
  
  return IsInside;
}

void BFLVorOperator::ResetVtxCache(void)
{ 
  fInVtxCache->Delete();
  fOutVtxCache->Delete();
}

Bool_t BFLVorOperator::VtxIsInsideNewPolyg(Int_t VtxID) const
{
// Returns TRUE if the vertex : VtxID is inside the Voronoi
// region that will be assigned to the new generator : igen.
// 
  Bool_t IsInside;

  // Check if Vtx is in our cache
  if (fInVtxCache->Exists(VtxID)) return kTRUE; 
  else if (fOutVtxCache->Exists(VtxID))  return kFALSE; 
  else { // not in cache.
    if (fVor->GetWeight(VtxID)) {  // If it is an ordinary point.

       TIntList * PolygsAround;  
       TObjArray * PointsOnCircle = new TObjArray(3); 

       // Get the number of polygons incident to vtx.
       // Should be 3, since we consistently avoid degeneracy.
       PolygsAround = RetrievePolygonsIncidentToVtx(VtxID);
       
       // Get the generators of these polygons.
       for(Int_t i=0; i < PolygsAround->NumberOfElements(); i++) { 
           Int_t GenID = fVor->PolygonToGenerator(PolygsAround->At(i));
           Float_t x = fVor->GetGeneratorX(GenID);
           Float_t y = fVor->GetGeneratorY(GenID);  
           PointsOnCircle->Add(new BFLNode(x,y));
       }

       // Check if the new generator is in the circle defined by
       // the points on the PointsOnCircle container.                       
       if (IsInTheCircle(PointsOnCircle,fCurrentGen) == 1) {
          IsInside = kTRUE;     
          fInVtxCache->AddLast(VtxID);
       } else {
          IsInside = kFALSE;
          fOutVtxCache->AddLast(VtxID);
       }
       
       // Explicitly free some allocated memory  
       delete PolygsAround;   
       PointsOnCircle->Delete();
       delete PointsOnCircle;

    } else { // point at infinity.     
       IsInside = kFALSE; 
       fOutVtxCache->AddLast(VtxID);
    } 
  } 
  
  return IsInside;  
} 


Int_t BFLVorOperator::IsInTheCircle(TObjArray * PointsOnCircle,
                                     BFLNode * point) const
{ 
// Checks the position of a point relative to the circle defined
// by the 3 points in the PointOnCircle container.
// It checks the determinant :
//
//             | 1   x1   y1   x1**2+y1**2 |
//             | 1   x2   y2   x2**2+y2**2 |
//        H =  | 1   x3   y3   x3**2+y3**2 |
//             | 1   x    y    x**2+y**2   |
//
// and returns:   1 if the point is within the circle
//                0 if the point is on the circle
//               -1 if the point is outside the circle.
//
// CAUTION: points (x1,y1), (x2,y2), (x3,y3) are placed in the xy 
// plane counterclockwisely ( when view from above this plane ).
//
  Int_t irow, InCircle;  
  Double_t det;
  BFLNode * PointOnCircle;
  
  TMatrix H = TMatrix(4,4);  

  // fill 1st column with 1's
  for(irow=0; irow<H.GetNrows(); irow++)  H(irow,0) = 1;
  
  TIter next(PointsOnCircle);
  irow = 0;
  while( (PointOnCircle = (BFLNode *)next())) {
     H(irow,1) = PointOnCircle->GetX();
     H(irow,2) = PointOnCircle->GetY();
//     H(irow,3) = TMath::Power(PointOnCircle->GetX(),2)+
//                         TMath::Power(PointOnCircle->GetY(),2);
     H(irow,3) = PointOnCircle->GetX()*PointOnCircle->GetX() +
                 PointOnCircle->GetY()*PointOnCircle->GetY();
     irow++;                         
  }
  
  // Fill the last row with point's distance from the origin
  H(3,1) = point->GetX();
  H(3,2) = point->GetY();
//  H(3,3) = TMath::Power(point->GetX(),2)+
//               TMath::Power(point->GetY(),2);
  H(3,3) = point->GetX()*point->GetX() + point->GetY()*point->GetY();
                       
  // Calculate the determinant(H)
  det = H.Determinant();                    
 
  if      (det == 0)  InCircle =  0;  // The point is on the circle
  else if (det  > 0)  InCircle = -1;  // The point is outside the circle
  else                InCircle =  1;  // The point is inside the circle
  
  // Explicitly free some allocated memory  
  delete PointOnCircle;

  return InCircle;
}

BFLNode BFLVorOperator::VtxPosition(TObjArray * PointsOnCircle) const
{
// Returns an BFLNode that holds the coordinates of the center of a
// circle that is described by the points on it ( PointsOnCirce 
// container of BFLNodes ).
// Actually, this is the way we calculate a vtx position as the center
// of the circle that passes through the generators of the three
// Voronoi regions around this vertex.
//
  Int_t i=0;
  Float_t A, B, C, Xvtx, Yvtx;
  Float_t x[3], y[3], r[3];
  BFLNode * PointOnCircle;
  
  TIter next(PointsOnCircle);
  while( (PointOnCircle = (BFLNode *)next())) {  
     x[i]=PointOnCircle->GetX();
     y[i]=PointOnCircle->GetY();
//     r[i]=TMath::Power(x[i],2)+TMath::Power(y[i],2);
     r[i]=x[i]*x[i] + y[i]*y[i];
     i++;
  }
  
  A = x[1]*y[2] - x[2]*y[1] - x[0]*y[2] +
      x[2]*y[0] + x[0]*y[1] - x[1]*y[0];
      
  B = y[1]*r[2] - y[2]*r[1] - y[0]*r[2] + 
      y[2]*r[0] + y[0]*r[1] - y[1]*r[0];
      
  C = - x[0]*r[1] - x[1]*r[2] - x[2]*r[0] + 
        x[0]*r[2] + x[1]*r[0] + x[2]*r[1];
      
  if(A !=0 ) {
     Xvtx = -B/(2*A);
     Yvtx = -C/(2*A);
     return BFLNode(Xvtx,Yvtx); 
  } else {
     return BFLNode(0.,0.);
  }     
}

BFLNode BFLVorOperator::FindNewVtx(Int_t EdgeID) const
{
// This method returns the new vertex coordinate on the edge: EdgeID.
// It defines the circle that passes through the generators of the
// left and right polygons ( with respect to EdgeID ) and the newly
// added generator, and then calls the VtxPosition method.
//  
  Int_t LeftPolygID, RightPolygID;
  Int_t LeftPolygGenID, RightPolygGenID;
  Float_t x,y;
  BFLNode vtx;
  TObjArray * PointsOnCircle = new TObjArray(3); 

  // Get the left and right polygon of the edge: EdgeID   
  LeftPolygID = fVor->GetLeftPolyg(EdgeID);
  RightPolygID = fVor->GetRightPolyg(EdgeID);
  
  // Get the IDs of the left and right polygon's generators
  LeftPolygGenID = fVor->GeneratorToPolygon(LeftPolygID);
  RightPolygGenID = fVor->GeneratorToPolygon(RightPolygID);
  
  // The above generators + the current generator define a circle
  // whose center is the new vertex on edge: EdgeID
  x = fVor->GetGeneratorX(LeftPolygGenID);
  y = fVor->GetGeneratorY(LeftPolygGenID);
  PointsOnCircle->Add(new BFLNode(x,y));

  x = fVor->GetGeneratorX(RightPolygGenID);
  y = fVor->GetGeneratorY(RightPolygGenID);
  PointsOnCircle->Add(new BFLNode(x,y));

  x = fCurrentGen->GetX();
  y = fCurrentGen->GetY();
  PointsOnCircle->Add(new BFLNode(x,y)); 

  vtx = VtxPosition(PointsOnCircle); // get the circle's center

  // Explicitly free some allocated memory
  PointsOnCircle->Delete(); 
  delete PointsOnCircle; 
    
  return vtx;
}

Float_t BFLVorOperator::GeneratorsDist(Int_t igen, Int_t jgen) const
{
// Returns the distance between the ith and jth generators.

   Float_t x = fVor->GetGeneratorX(igen) - fVor->GetGeneratorX(jgen);
   Float_t y = fVor->GetGeneratorY(igen) - fVor->GetGeneratorY(jgen);
  return TMath::Sqrt(
//         TMath::Power(fVor->GetGeneratorX(igen)- 
//                      fVor->GetGeneratorX(jgen),2) +
//         TMath::Power(fVor->GetGeneratorY(igen)- 
//                      fVor->GetGeneratorY(jgen),2) );  
     x*x + y*y);
}

Float_t BFLVorOperator::DistanceFrom(Int_t igen) const
{
// Returns the distance between the ith and jth generators.
// 

   Float_t x = fVor->GetGeneratorX(igen) - fCurrentGen->GetX();
   Float_t y = fVor->GetGeneratorY(igen) - fCurrentGen->GetY();
  return TMath::Sqrt(
//       TMath::Power(fVor->GetGeneratorX(igen)- fCurrentGen->GetX(),2) +
//       TMath::Power(fVor->GetGeneratorY(igen)- fCurrentGen->GetY(),2) );  
     x*x + y*y);
}

Int_t BFLVorOperator::MoveToNextPolygon(Int_t PolygID, Int_t EdgeID) const
{
// Returns the polygon that shares the Voronoi Edge: EdgeID  other 
// than the Voronoi Polygon: PolygID.

  Int_t LeftPolygID, RightPolygID, NextPolygID;
  
  LeftPolygID = fVor->GetLeftPolyg(EdgeID);
  RightPolygID = fVor->GetRightPolyg(EdgeID);

  if(LeftPolygID != PolygID) NextPolygID = LeftPolygID;
  else NextPolygID = RightPolygID;

  return NextPolygID;
}

Int_t BFLVorOperator::FindCurrentPolygon(Int_t GuessedPolygID) const
{
// Returns the Voronoi polygon that contains the added generator.
// In a straightforward implementation the computational complexitity
// is ~ N, but within a ~ N loop it gives a computational complexity
// of ~ N**2. This is what we want to avoid since N ~ 100k.
// 
// igenerator is the id of the added generator 
        
  Int_t i, NearestGenID, Nedges, AdjacentPolygonID;
  Int_t PolygID, EdgeID, AdjGen;
  Bool_t IsNearest;
  Float_t Dist, DistMin;
  TIntList * SurroundingEdges;

  // Get a guess of the nearest generator (if any). 
  if(GuessedPolygID != -1 ) NearestGenID = GuessedPolygID;
  else NearestGenID = 1;

  // Distance between current and "nearest" generators.
  DistMin = DistanceFrom(NearestGenID);
  
  // Look for a closer generator than our guess. 
  do {  
    IsNearest = kTRUE;

    // Voronoi polygon corresponding to "nearest" generator.    
    PolygID = fVor->GeneratorToPolygon(NearestGenID);

    // Get the edges surrounding the polygon that belongs to the
    // closest generator.    
    //cout << "...surrounding edges" << endl;
    SurroundingEdges = RetrieveEdgesSurrPolygon(PolygID);
   
    // Loop over surrounding edges.    
    //cout << "... loop over surrounding edges" << endl;
    Nedges = SurroundingEdges -> NumberOfElements();
    for(i=0; i<Nedges; i++) {

       // Get the edge and the adjacent polygon (other than PolygID).
       EdgeID = SurroundingEdges -> At(i); 
       AdjacentPolygonID = MoveToNextPolygon(PolygID,EdgeID); 
    
       if(AdjacentPolygonID !=0 ) {
           // Get the generator of the adjacent polygon
           // and its coordinates.
           AdjGen = fVor->PolygonToGenerator(AdjacentPolygonID);

           // Distance between the current and adjacent polygon's 
           // generators.
           //cout << "Adjacent Gen = " << AdjGen << endl;
           Dist = DistanceFrom(AdjGen);
           
           // Decide if this is a closer generator.
           if ( Dist < DistMin - 0.0001) {
                DistMin = Dist;
                NearestGenID = AdjGen;
                IsNearest = kFALSE;
           }       
       } // if polygid != 0       
    } // from loop over edges
    
    // Explicitly free some allocated memory
    delete SurroundingEdges;
   
  } while(!IsNearest);
  
  return fVor->GeneratorToPolygon(NearestGenID);
}

---