Back to index

PHGeometron.C

 
//--------------------------------------------------------------------------------------- 
//  
// Includes:  PHGeometron.h 
// 
// Created by:  Jane M. Burward-Hoy and Federica Messer 
// 
// Purpose:  singleton for Geometry transformations 
//  
// Last update:  10/19/99 
// 
//--------------------------------------------------------------------------------- 
#include "PHGeometron.h" 
#include <stdlib.h> 
#include "phool.h" 
 
PHGeometron* PHGeometron::_instance = 0; 
 
PHGeometron::PHGeometron()  
{ 
   
} 
 
PHGeometron::~PHGeometron()  
{ 
 
} 
 
PHGeometron * PHGeometron::instance() 
{ 
   if (!_instance) 
      _instance = new PHGeometron; 
   return _instance; 
} 
 
 
void PHGeometron::cartesianToCylindrical(const PHPoint &cart, PHCylPoint& cyl) 
{ 
  double r; 
  PHAngle phi; 
 
  r = sqrt(cart.getX()*cart.getX() + cart.getY()*cart.getY()); 
   
  if (cart.getX() != 0 ) { 
    phi = atan2(cart.getY(), cart.getX()); 
  }else { 
    if (cart.getY() > 0.0) { 
      phi = Pi/2.; 
    }else if (cart.getY() < 0.0) { 
      phi = -Pi/2.; 
    }else { 
      phi = 0; 
    } 
  } 
   
  cyl.setR(r); 
  cyl.setPhi(phi); 
  cyl.setZ(cart.getZ()); 
   
  cout << "check sign and values of phi " << endl; 
   
} 
 
void PHGeometron::cylindricalToCartesian(const PHCylPoint& cyl, PHPoint& cart) 
{ 
  double x,y; 
 
  x = cyl.getR()*cos(cyl.getPhi()); 
  y = cyl.getR()*sin(cyl.getPhi()); 
   
  cart.setX(x); 
  cart.setY(y); 
  cart.setZ(cyl.getZ()); 
  
  cout << "check sign and values with Jeff " << endl; 
   
} 
 
void PHGeometron::cartesianToSpherical(const PHPoint &cart, PHSphPoint& sph) 
{ 
  double r; 
  PHAngle phi,theta; 
 
  r =  sqrt(cart.getX()*cart.getX() + cart.getY()*cart.getY() + cart.getZ()*cart.getZ()); 
  if (cart.getX() != 0) { 
    phi = atan2(cart.getY(), cart.getX()); 
  } else { 
    if (cart.getY() > 0 ) { 
      phi = Pi/2; 
    }else if (cart.getY() < 0) { 
      phi = -Pi/2; 
    }else { 
      phi = 0; 
    } 
  } 
 
  if ( r != 0 ) { 
    theta = acos(cart.getZ()/r); 
  }else { 
    theta = 0; 
  } 
 
  sph.setR(r); 
  sph.setPhi(phi); 
  sph.setTheta(theta); 
  cout << "check sign and values with Jeff " << endl; 
  cout << " NB: jeff invert theta and Phi" << endl; 
 
} 
 
void PHGeometron::sphericalToCartesian(const PHSphPoint& sph, PHPoint& cart) 
{ 
  cart.setX(sph.getR()*sin(sph.getTheta())*cos(sph.getPhi())); 
  cart.setY(sph.getR()*sin(sph.getTheta())*sin(sph.getPhi())); 
  cart.setZ(sph.getR()*cos(sph.getTheta())); 
} 
 
void PHGeometron::cylindricalToSpherical(const PHCylPoint& cyl,PHSphPoint &sph) 
{ 
 PHPoint cart; 
 cylindricalToCartesian(cyl,cart); 
 cartesianToSpherical(cart,sph); 
   
} 
 
void PHGeometron::sphericalToCylindrical(const PHSphPoint &sph,PHCylPoint& cyl) 
{ 
  PHPoint cart; 
  sphericalToCartesian(sph,cart); 
  cartesianToCylindrical(cart,cyl); 
   
} 
 
double PHGeometron::distancePointToPoint(const PHPoint& c1, const PHPoint &c2) 
{ 
 return sqrt( pow(c1.getX()-c2.getX(),2) + 
	      pow(c1.getY()-c2.getY(),2) + 
	      pow(c1.getZ()-c2.getZ(),2) );   
} 
 
double PHGeometron::distancePointToPoint2D(const PHPoint& c1, const PHPoint &c2) 
{ 
  PHPoint c1point2D = c1;    
  PHPoint c2point2D = c2; 
  c1point2D.setZ(0.); 
  c2point2D.setZ(0.); 
  return distancePointToPoint(c1point2D,c2point2D); 
} 
 
double PHGeometron::distancePointToPoint(const PHSphPoint& s1, const PHSphPoint &s2) 
{ 
 PHPoint c1 = s1; 
 PHPoint c2 = s2; 
 return distancePointToPoint(c1,c2); 
} 
 
double PHGeometron::distancePointToPoint(const PHCylPoint& cy1, const PHCylPoint &cy2) 
{ 
 PHPoint c1 = cy1; 
 PHPoint c2 = cy2; 
 return  distancePointToPoint(c1,c2); 
} 
 
     
double  PHGeometron::dot(const PHVector& v1, const PHVector& v2) 
{ 
  double product; 
  product = v1.getX()*v2.getX() + v1.getY()*v2.getY() + v1.getZ()*v2.getZ(); 
  return product; 
} 
 
PHVector PHGeometron::cross(const PHVector& v1, const PHVector& v2) 
{ 
  
  double x, y, z; 
 
  x = (v1.getY()*v2.getZ() - v1.getZ()*v2.getY()); 
  y = (v1.getZ()*v2.getX() - v1.getX()*v2.getZ()); 
  z = (v1.getX()*v2.getY() - v1.getY()*v2.getX()); 
 
  PHVector out(x, y, z); 
  return out; 
} 
 
PHAngle PHGeometron::angle(const PHVector& v1, const PHVector& v2) 
{ 
  PHAngle phi; 
  if ((v1.length() != 0.0) && (v2.length()!=0.0)) { 
    phi = acos( v1.dot(v2)/sqrt(v1.lengthSqr()* v2.lengthSqr())); 
  }else { 
    cout << "PHGeometron::angle -- Length of one vector is zero "<< endl; 
  } 
   
  return phi;   
} 
 
PHBoolean PHGeometron::intersectionLinePlane(const PHLine& line, const PHPlane& plane, PHPoint& cross) 
{ 
  PHVector lineVector  = line.getDirection(); 
  PHVector planeVector = plane.getNormal(); 
  PHPoint  linePoint   = line.getBasepoint(); 
  PHPoint  planePoint  = plane.getOrigin(); 
 
  PHVector baseVector = planePoint - linePoint; // vector connecting the basepoints 
 
  double distance,denominator; 
   
  denominator = planeVector.dot(lineVector); 
 
  if ( denominator != 0) { // if line and plane are not parallel 
    distance = (planeVector.dot(baseVector))/(planeVector.dot(lineVector)); 
    cross = linePoint + (PHPoint)(lineVector*distance); 
    return True; 
  }else {                   // if line and plane are parallel 
    return False; 
  } 
} 
 
PHBoolean PHGeometron::intersectionLinePlane2D(const PHLine& line, const PHPlane& plane, PHPoint& cross) 
{ 
  PHVector lineVector  = line.getDirection(); 
  PHVector planeVector = plane.getNormal(); 
  PHPoint  linePoint   = line.getBasepoint(); 
  PHPoint  planePoint  = plane.getOrigin(); 
 
  lineVector.setZ(0.); 
  planeVector.setZ(0.); 
  linePoint.setZ(0.); 
  planePoint.setZ(0.); 
 
  PHPlane plane2D(planePoint, planeVector); 
  PHLine  line2D(linePoint, lineVector); 
 
  if (intersectionLinePlane(line2D, plane2D,cross)) { 
    return True; 
  }else { 
    return False; 
  } 
} 
/* 
PHBoolean PHGeometron::intersectionLineLine2D(const PHLine& line1, const PHLine &line2, PHPoint& cross) 
{ 
  PHVector vector1 = line1.getDirection(); 
  PHVector vector2 = line2.getDirection(); 
  PHPoint  point1  = line1.getBasepoint(); 
  PHPoint  point2  = line2.getBasepoint(); 
 
  double m1, m2;                                 //slopes of lines 
  double b1, b2;                                //y-intercepts of lines 
  double denominator1, denominator2;      
 
  vector1.setZ(0.); 
  vector2.setZ(0.); 
  point1.setZ(0.); 
  point2.setZ(0.); 
 
  if ( (vector1.length() == 0.0) ||  (vector2.length() == 0.0) ) { 
	  cout << "ERROR: Null line in PHGeometron::intersectionLineLine2D"  
	       << endl; 
	  return False; } //Null line 
  else { 
    denominator1 = vector1.getX(); 
    denominator2 = vector2.getX(); 
    if ( (denominator1 != 0.0) && (denominator2 != 0.0) ) { 
      m1 = vector1.getY()/vector1.getX(); 
      b1 = -m1*point1.getX() + point1.getY(); 
      m2 = vector2.getY()/vector2.getX(); 
      b2 = -m2*point2.getX() + point2.getY(); 
      if ( (m1 != 0.0) && (m2 !=0.0) ) 
	if ( ( m1 != m2 ) && ( m1 != -m2 ) ) { 
	  cout << "L1 and L2 are not parallel" << endl; 
	  cross.setX((b2 - b1)/(m1 - m2)); 
	  cross.setY( (m1*b2 - m2*b1)/(m1 - m2)); 
	  return True;  } //Not parallel  
	else {   //Parallel lines, m1 == m2 or m1 = -m2   
	  cout << "L1 and L2 are parallel" << endl; 
	  return False; }  
      else {    //Horizontal case 
	if ( ( m1 == 0.0 ) && (m2 != 0.0) ) {              //vector1 is horizontal 
	    cout << "L1 is Horizontal" << endl; 
	    cross.setX( (point2.getY() - b2) / m2 ); 
	    cross.setY( point1.getY() ); 
	    return True; }                                
	if ( ( m2 == 0.0) && ( m1 != 0.0) ) {             //vector2 is horizontal 
	    cout << "L2 is Horizontal" << endl; 
	    cross.setX( (point1.getY() - b1) /m1 ); 
	    cross.setY( point2.getY() ); 
	    return True;}  
        if ( (m1 == 0.0) && (m2 == 0.0) ) { 
	    cout << "L1 and L2 are parallel" << endl; 
	    return False;} 
    } }  // end of no zero denominator case 
    else         //Vertical line case 
      if (denominator1 == 0.0) {   //vector1 is vertical 
	cout << "L1 is Vertical" << endl; 
	m2 = vector2.getY()/vector2.getX(); 
	b2 = -m2*point2.getX() + point2.getY(); 
	cross.setX(point1.getX()); 
	cross.setY(m2*point1.getX() + b2); 
	return True; } 
      else  {                      //vector2 is vertical 
	cout << "L2 is Vertical" << endl; 
	m1 = vector1.getY()/vector1.getX(); 
	b1 = -m1*point1.getX() + point1.getY(); 
	cross.setX(point2.getX()); 
	cross.setY(m1*point2.getX() + b1); 
	return True; } 
  }  //the main "else" closure 
} 
*/ 
 
 
/* 
 
PHBoolean PHGeometron::intersectionLineLineOnPlane(const PHLine& line1, const PHLine& line2, const PHPlane& plane, PHPoint& cross) 
{ 
  
  PHLine projected1, projected2; 
 
  projectLinePlane(line1,plane,projected1); 
  projectLinePlane(line2,plane,projected2); 
 
  PHVector vector1 = projected1.getDirection(); 
  PHVector vector2 = projected2.getDirection(); 
 
  PHPoint point1 = projected1.getBasepoint(); 
  PHPoint point2 = projected2.getBasepoint(); 
 
  
  //The code that follows is just like intersectionLineLine2D( ... ) above 
  //Note: we are now in the plane's coordinate system.  The origin is the 
  //plane's origin. 
 
  double m1, m2;                                 //slopes of lines 
  double b1, b2;                                //y-intercepts of lines 
  double denominator1, denominator2;      
 
 
  if ( (vector1.length() == 0.0) ||  (vector2.length() == 0.0) ) { 
	  cout << "ERROR: Null line in PHGeometron::intersectionLineLine2D"  
	       << endl; 
	  return False; } //Null line 
  else { 
    denominator1 = vector1.getX(); 
    denominator2 = vector2.getX(); 
    if ( (denominator1 != 0.0) && (denominator2 != 0.0) ) { 
      m1 = vector1.getY()/vector1.getX(); 
      b1 = -m1*point1.getX() + point1.getY(); 
      m2 = vector2.getY()/vector2.getX(); 
      b2 = -m2*point2.getX() + point2.getY(); 
      if ( (m1 != 0.0) && (m2 !=0.0) ) 
	if ( ( m1 != m2 ) && ( m1 != -m2 ) ) { 
	  cout << "L1 and L2 are not parallel" << endl; 
	  cross.setX((b2 - b1)/(m1 - m2)); 
	  cross.setY( (m1*b2 - m2*b1)/(m1 - m2)); 
	  return True;  } //Not parallel  
	else {   //Parallel lines, m1 == m2 or m1 = -m2   
	  cout << "L1 and L2 are parallel" << endl; 
	  return False; }  
      else    //Horizontal case 
	if ( m1 == 0.0 ) {              //vector1 is horizontal 
	    cout << "L1 is Horizontal" << endl; 
	    cross.setX( (point2.getY() - b2) / m2 ); 
	    cross.setY( point1.getY() ); 
	    return True; } 
	else {                  //vector2 is horizontal 
	    cout << "L2 is Horizontal" << endl; 
	    cross.setX( (point1.getY() - b1) /m1 ); 
	    cross.setY( point2.getY() ); 
	    return True;}    } // end of no zero denominator case 
    else         //Vertical line case 
      if (denominator1 == 0.0) {   //vector1 is vertical 
	cout << "L1 is Vertical" << endl; 
	m2 = vector2.getY()/vector2.getX(); 
	b2 = -m2*point2.getX() + point2.getY(); 
	cross.setX(point1.getX()); 
	cross.setY(m2*point1.getX() + b2); 
	return True; } 
      else  {                      //vector2 is vertical 
	cout << "L2 is Vertical" << endl; 
	m1 = vector1.getY()/vector1.getX(); 
	b1 = -m1*point1.getX() + point1.getY(); 
	cross.setX(point2.getX()); 
	cross.setY(m1*point2.getX() + b1); 
	return True; } 
  }  //the main "else" closure 
 
 
} 
 
*/ 
double PHGeometron::distanceLinePoint(const PHLine& line, const PHPoint &point) 
{ 
  PHPlane myplane(point, line.getDirection()); // plane perp. to line 
  PHPoint cross; 
  double distance; 
   
  if (intersectionLinePlane(line, myplane, cross)) { 
     distance = (PHVector(point - cross)).length(); 
     return distance; 
  } else { 
    cout << "ERROR in PHGeometron::distanceLinePoint() " << endl; 
    // if an error, plane and line are  parallel or vector length = 0. 
    // -->check into intersectionLinePlane 
    return 0.0; 
  } 
} 
 
double PHGeometron::distanceLinePoint2D(const PHLine& line, const PHPoint &point) 
{ 
  PHPoint point2D = point; 
  point2D.setZ(0.); 
 
  PHLine line2D = line; 
  line.getDirection().setZ(0.); 
  line.getBasepoint().setZ(0.); 
   
  return distanceLinePoint(line2D, point2D); 
   
} 
 
 
double PHGeometron::distanceLineLine(const PHLine& line1, const PHLine &line2) 
{ 
  double distance = 0.; 
  PHPoint cross; 
  PHVector vector1 = line1.getDirection(); 
  PHVector vector2 = line2.getDirection(); 
  PHAngle angleBetweenLines = vector1.angle(vector2); 
  
  if ((double)angleBetweenLines == 0.) { // if lines are not parallel 
    PHPlane myplane(line2.getBasepoint(),line2.getDirection()); 
    if (intersectionLinePlane(line1, myplane,cross)) { 
      distance = (PHVector(cross - line2.getBasepoint())).length(); 
    }else { 
      cout <<  "ERROR in PHGeometron:: distanceLineLine" << endl; 
    } 
  } else { 
     PHVector connection = vector1.cross(vector2); 
     PHVector normal = vector1.cross(connection); 
     PHPlane myplane(line1.getBasepoint(),normal); 
     if( intersectionLinePlane(line2,myplane,cross)) { 
       distance = distanceLinePoint(line1,cross); 
     }else{ 
      cout <<  "ERROR in PHGeometron:: distanceLineLine" << endl; 
     } 
  } 
  return distance; 
} 
 
 
short PHGeometron::intersectionLineSphere(const PHLine& line, const PHSphere& sph, 
					      PHPoint& cross1, PHPoint& cross2) 
{ 
 
  PHPoint  linePoint  = line.getBasepoint(); 
  PHVector lineVector = line.getDirection(); 
  lineVector.normalize(); 
 
  PHPoint center = sph.getCenter(); 
  double radius  = sph.getRadius(); 
 
  // Calculate the coefficient of the quadratic equation 
  double a = lineVector.dot(lineVector); 
  double b = 2.0 *(lineVector.dot((linePoint-center))); 
  PHVector tmpVector = linePoint - center; 
  double c = tmpVector.dot(tmpVector) - radius*radius; 
  double d = b*b - 4.0*a*c; 
 
  cout << " d ;" << d << endl; 
  // if d <0 or a = 0, there is no intercept 
  if (d <0.0 || a==0) return 0; 
 
  // if d = 0; there is only one intercept at  -b/2a 
  if ( d == 0) { 
    cross1 = linePoint + PHPoint(lineVector*(-b/(2*a))); 
    return 1; 
  } 
  // if d > 0, there are 2 intercept 
  if (d > 0) { 
    cross1 = linePoint +PHPoint(lineVector*((-b + sqrt(d))/(2.*a))); 
    cross2 = linePoint +PHPoint(lineVector*((-b - sqrt(d))/(2.*a))); 
    return 2; 
  } 
   
  return True; 
} 
 
short PHGeometron::intersectionLineCylinder(const PHLine& line, const PHCylinder& cyl, 
						PHPoint& cross1, PHPoint& cross2) 
{ 
  double uC,uL; 
  double a ,b,c,d ; 
  double u1,u2; 
  PHPoint pCyl; 
  PHPoint  linePoint  = line.getBasepoint(); 
  PHVector lineVector = line.getDirection(); 
  lineVector.normalize(); 
 
   
  PHPoint center = cyl.getCenter(); 
  PHVector axis  = cyl.getAxis(); 
  double radius = cyl.getRadius(); 
   
  // special case  
  if(lineVector.getZ() == 0.) { 
    if(axis.getZ() !=0) { 
      uC = (linePoint.getZ() - center.getZ())/axis.getZ(); 
      pCyl = center + (PHPoint)(axis*uC); 
      a = lineVector.dot(lineVector); 
      b = 2.0*(lineVector.dot((PHVector)(linePoint - pCyl))); 
      c = (PHVector(linePoint-pCyl)).dot((PHVector(linePoint-pCyl))) - radius*radius 
	  - (linePoint.getZ() - pCyl.getZ())*(linePoint.getZ() - pCyl.getZ()); 
      d = b*b -4.*a*c; 
 
      if (a == 0 || d < 0) return 0;  
      u1 = (-b + sqrt(d))/(2.*a); 
      u2 = (-b - sqrt(d))/(2.*a); 
			   
    }else {  // there is no intercept 
      return 0;  
    } 
  } else { // lineVector.getZ()  
    double a1,a2,a3,a4,a5; 
    a1 = (linePoint.getX() - center.getX())*(linePoint.getX() - center.getX()) + 
         (linePoint.getY() - center.getY())*(linePoint.getY() - center.getY()) -radius*radius; 
    a2 = 2.0*(lineVector.getX()*(linePoint.getX()-center.getX()) + 
	      lineVector.getY()*(linePoint.getY()-center.getY())); 
    a3 = 2.0*(axis.getX()*(center.getX() -linePoint.getX()) + 
	      axis.getY()*(center.getY() -linePoint.getY())); 
    a4 = lineVector.getX()*lineVector.getX() + lineVector.getY()*lineVector.getY(); 
    a5 = axis.getX()*axis.getX() +  axis.getY()*axis.getY(); 
 
    a = a4*(axis.getZ()/lineVector.getZ())*(axis.getZ()/lineVector.getZ()) + a5; 
    b = a2*(axis.getZ()/lineVector.getZ()) +2.0*a4*(axis.getZ()/lineVector.getZ()) +a3; 
    c = (a4*(center.getZ()/linePoint.getZ())*(center.getZ()/linePoint.getZ()))/(lineVector.getZ()*lineVector.getZ()) + 
        (a2*(center.getZ()/linePoint.getZ()))/(lineVector.getZ()) +a1; 
    
    d = b*b -4.0*a*c; 
   
    if (a == 0 || d < 0) return 0; 
 
    double uc1,uc2; 
    uc1 = (-b +sqrt(d))/(2.0*a); 
    uc2 = (-b -sqrt(d))/(2.0*a); 
    u1 = (center.getZ() - linePoint.getZ() + uc1*axis.getZ())/lineVector.getZ(); 
    u2 = (center.getZ() - linePoint.getZ() + uc2*axis.getZ())/lineVector.getZ(); 
  } 
 
 
  // if d = 0 , there is one intersection 
 
  if (d == 0) { 
    cross1 = linePoint + (PHPoint)lineVector*(u1); 
    cross2 = linePoint + (PHPoint)lineVector*(u2); 
    return 1; 
  } 
  // if d > 0 , there are 2 intersection 
  if (d > 0) { 
    cross1 = linePoint + (PHPoint)lineVector*(u1); 
    cross2 = linePoint + (PHPoint)lineVector*(u2); 
    return 2; 
  } 
   
   return True; 
} 
 
/* 
 
short  PHGeometron::lineLineIntersection(const PHLine&l1, const PHLine &l2, PHPoint &) 
{ 
} 
*/ 
 
 
PHLine PHGeometron::projectLineIntoPlane(const PHLine& line, const PHPlane& plane) 
{ 
  PHVector vector = line.getDirection(); 
  PHPoint  point  = line.getBasepoint(); 
  PHVector normal = plane.getNormal(); 
  PHPoint origin  = plane.getOrigin(); 
   
  // Project direction vector of line onto plane 
   
  PHVector v_projected = vector - normal*(vector.dot(normal)); 
 
  //Project basepoint of line onto plane 
 
 
  PHPoint p_projected, p_ontonormal; 
  PHVector newvector = point - origin; 
  
  p_ontonormal = normal*newvector.dot(normal); 
  p_projected = point - p_ontonormal; 
 
  PHLine newline(p_projected,v_projected); 
   
  return newline;  
} 
 
PHPoint PHGeometron::transformPoint(const PHFrame &F1, const PHPoint &p1, const PHFrame &F2) 
  // Transform a point in one coordinate frame, F2 to another coordinate frame, F2. 
  // Return the transformed point. 
{ 
  //Translate the point: 
  PHPoint origin_distance = F2.getOrigin() - F1.getOrigin(); 
  PHVector translate; 
  translate = p1 - origin_distance; 
 
  //Calculate the direction cosines: 
  PHVector dir_cosX; 
  PHVector dir_cosY; 
  PHVector dir_cosZ; 
 
  PHVector u1, v1, w1; 
  PHVector u2, v2, w2; 
 
  u1 = F1.getU();  v1 = F1.getV(); w1 = F1.getW(); 
  u2 = F2.getU();  v2 = F2.getV(); w2 = F2.getW(); 
 
  u1.normalize(); v1.normalize(); w1.normalize(); 
  u2.normalize(); v2.normalize(); w2.normalize(); 
 
  dir_cosX.setX( u2.dot( u1 ) ); 
  dir_cosX.setY( u2.dot( v1 ) ); 
  dir_cosX.setZ( u2.dot( w1 ) ); 
   
  dir_cosY.setX( v2.dot( u1 ) ); 
  dir_cosY.setY( v2.dot( v1 ) ); 
  dir_cosY.setZ( v2.dot( w1 ) ); 
  
  dir_cosZ.setX( w2.dot( u1 ) ); 
  dir_cosZ.setY( w2.dot( v1 ) ); 
  dir_cosZ.setZ( w2.dot( w1 ) ); 
    
  //Rotate the translated point: 
  PHPoint rotate; 
 
  rotate.setX( translate.dot(dir_cosX) ); 
  rotate.setY( translate.dot(dir_cosY) ); 
  rotate.setZ( translate.dot(dir_cosZ) ); 
 
  return rotate; 
 
} 
 
PHVector PHGeometron::transformVector(const PHFrame &F1, const PHVector &v, const PHFrame &F2) 
{ 
  PHPoint new_cast_p; 
  PHPoint cast_p; 
  PHVector transformed; 
  
  cast_p = v; 
  new_cast_p = transformPoint(F1, cast_p, F2); 
   
  transformed = new_cast_p; 
  return transformed; 
} 
 
PHLine PHGeometron::transformLine(const PHFrame &F1, const PHLine &line, const PHFrame &F2) 
{ 
  PHPoint  point  = line.getBasepoint(); 
  PHVector vector = line.getDirection(); 
 
  PHPoint  transformed_p = transformPoint(F1, point, F2); 
  PHVector transformed_v = transformVector(F1, vector, F2); 
 
  PHLine transformed(transformed_p,transformed_v); 
 
  return transformed; 
} 
 
 
// PHBoolean PHGeometron::intersectionLineLineOnPlane(const PHLine& line1, const PHLine& line2, const PHPlane& plane, PHPoint& cross) 
// { 
//   PHLine projected1 = projectLinePlane(line1,plane); 
//   PHLine projected2 = projectLinePlane(line2,plane); 
 
//   PHVector vector1 = projected1.getDirection(); 
//   PHVector vector2 = projected2.getDirection(); 
 
//   PHPoint point1 = projected1.getBasepoint(); 
//   PHPoint point2 = projected2.getBasepoint(); 
 
  
//   //The code that follows is just like intersectionLineLine2D( ... ) above 
//   //Note: we are now in the plane's coordinate system.  The origin is the 
//   //plane's origin. 
 
//   double m1, m2;                                 //slopes of lines 
//   double b1, b2;                                //y-intercepts of lines 
//   double denominator1, denominator2;      
 
 
//   if ( (vector1.length() == 0.0) ||  (vector2.length() == 0.0) ) { 
// 	  cout << "ERROR: Null line in PHGeometron::intersectionLineLine2D"  
// 	       << endl; 
// 	  return False; } //Null line 
//   else { 
//     denominator1 = vector1.getX(); 
//     denominator2 = vector2.getX(); 
//     if ( (denominator1 != 0.0) && (denominator2 != 0.0) ) { 
//       m1 = vector1.getY()/vector1.getX(); 
//       b1 = -m1*point1.getX() + point1.getY(); 
//       m2 = vector2.getY()/vector2.getX(); 
//       b2 = -m2*point2.getX() + point2.getY(); 
//       if ( (m1 != 0.0) && (m2 !=0.0) ) 
// 	if ( ( m1 != m2 ) && ( m1 != -m2 ) ) { 
// 	  cout << "L1 and L2 are not parallel" << endl; 
// 	  cross.setX((b2 - b1)/(m1 - m2)); 
// 	  cross.setY( (m1*b2 - m2*b1)/(m1 - m2)); 
// 	  return True;  } //Not parallel  
// 	else {   //Parallel lines, m1 == m2 or m1 = -m2   
// 	  cout << "L1 and L2 are parallel" << endl; 
// 	  return False; }  
//       else    //Horizontal case 
// 	if ( m1 == 0.0 ) {              //vector1 is horizontal 
// 	    cout << "L1 is Horizontal" << endl; 
// 	    cross.setX( (point2.getY() - b2) / m2 ); 
// 	    cross.setY( point1.getY() ); 
// 	    return True; } 
// 	else {                  //vector2 is horizontal 
// 	    cout << "L2 is Horizontal" << endl; 
// 	    cross.setX( (point1.getY() - b1) /m1 ); 
// 	    cross.setY( point2.getY() ); 
// 	    return True;}    } // end of no zero denominator case 
//     else         //Vertical line case 
//       if (denominator1 == 0.0) {   //vector1 is vertical 
// 	cout << "L1 is Vertical" << endl; 
// 	m2 = vector2.getY()/vector2.getX(); 
// 	b2 = -m2*point2.getX() + point2.getY(); 
// 	cross.setX(point1.getX()); 
// 	cross.setY(m2*point1.getX() + b2); 
// 	return True; } 
//       else  {                      //vector2 is vertical 
// 	cout << "L2 is Vertical" << endl; 
// 	m1 = vector1.getY()/vector1.getX(); 
// 	b1 = -m1*point1.getX() + point1.getY(); 
// 	cross.setX(point2.getX()); 
// 	cross.setY(m1*point2.getX() + b1); 
// 	return True; } 
//   }  //the main "else" closure 
// } 
 
PHLine PHGeometron::translateLinePoint(const PHLine &line, const PHPoint &p) 
{ 
  // Translate the line to the point p.  Simply set the basepoint of the translated 
  // line equal to p, but keep the direction of the original line. 
 
 PHLine translatedLine; 
 
 translatedLine.setBasepoint(p); 
 translatedLine.setDirection(line.getDirection()); 
 
 return translatedLine; 
 
} 
 
PHLine PHGeometron::translateLineVector(const PHLine &line, const PHVector &vector) 
{ 
 
  // Move the line (without changing its direction) to the point where the vector is pointing (the tip). 
 
  PHLine translatedLine; 
 
  translatedLine.setBasepoint( (PHPoint)vector + line.getBasepoint() ); 
  translatedLine.setDirection( line.getDirection() ); 
 
  return translatedLine; 
 
} 
 
double PHGeometron::lengthPolyLine(const PHPolyLine &PolyLine) 
{ 
  
  double totalLength = 0.0; 
  size_t first = 0; 
  size_t last = PolyLine.numPoints() - 1; 
  //NEEDS TO BE DEBUGGED! 
  // totalLength = PolyLine.PHPolyLine::length(first,last); 
 
  return (double)last; 
 
} 
 
 
// PHBoolean PHGeometron::intersectionPolyLinePlane(const PHPolyLine &PolyLine, const PHPlane &Plane, PHPoint &Point) 
// { 
//  PHBoolean ifIntersect = False; 
//  PHLine segment; 
//  PHPoint potentialPoint; 
 
//  for (size_t i = 0; i < PolyLine.numPoints(); i++) 
//    { 
 
//      PHPoint b( PolyLine.getPoint(i)->getX(), PolyLine.getPoint(i)->getY(), PolyLine.getPoint(i)->getZ()); 
 
//      PHVector v( (PolyLine.getPoint(i+1)->getX() - PolyLine.getPoint(i)->getX()), 
// 	        (PolyLine.getPoint(i+1)->getY() - PolyLine.getPoint(i)->getY()), 
// 		(PolyLine.getPoint(i+1)->getZ() - PolyLine.getPoint(i)->getZ()) ); 
      
//      segment.setBasepoint(b); 
//      segment.setDirection(v); 
 
//      if (intersectionLinePlane(segment,Plane,potentialPoint) ) 
 
//        if ( distancePointToPoint(potentialPoint,segment.getBasepoint()) < segment.length() ) 
// 	 { 
// 	   Point = potentialPoint; 
// 	   ifIntersect = True; 
// 	 } 
//    } 
//  return ifIntersect; 
// }  
 
 
 
 
 
 

Back to index