geometry.cpp

//
//  file    : geometry.cpp
//  author  : jh
//  date    : 15 jul 2002
//  changed : 25 oct 2002
//


#include <stdio.h>

#include "geometry.h"
#include "Vertex.h"


double
geometry::area(Vertex* v1, Vertex* v2, Vertex* v3)
{

  double a = (0.5 * ((v2->x() * v3->y()  -  v3->x() * v2->y()) +
                     (v3->x() * v1->y()  -  v1->x() * v3->y()) +    
                     (v1->x() * v2->y()  -  v2->x() * v1->y())) );  

  return a;

}



double 
geometry::distance(const point& p1, const point& p2)
{
  double d = sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) );
  return d;
} 
 


double
geometry::distance(const line& l, const point& p)
{
  double bX = l.p2.x - l.p1.x;
  double bY = l.p2.y - l.p1.y;
    
  double pX = p.x - l.p1.x;
  double pY = p.y - l.p1.y;
    
  double bxp = bX * pY - bY * pX;
    
  return fabs(bxp/distance(l.p1,l.p2));
  
}


bool
geometry::intersect(const line& l1, const line& l2)
{
  double CEPS = 1e-5;

#ifdef _DEBUG_
  fprintf(stderr, "<geometry::intersect()>\n");
  fprintf(stderr, " check (%f/%f)-(%f/%f) x \n",
          l1.p1.x, l1.p1.y,  l1.p2.x, l1.p2.y );
  fprintf(stderr, "       (%f/%f)-(%f/%f) \n",
          l2.p1.x, l2.p1.y,  l2.p2.x, l2.p2.y );
#endif
  
  double xmin, xmax, ymin, ymax;
  // double x1,x2,y1,y2;     
  double x21,x31,x41,x43; 
  double y21,y31,y41,y43; 
  double d0,d1,d2;        
  double s,t;         


  if( distance(l1.p1, l1.p2) < CEPS ) {
#ifdef _DEBUG_
    fprintf(stderr, "(1) return <false>\n");
#endif
    return false;
  }


  if( distance(l2.p1, l2.p2) < CEPS ) {
#ifdef _DEBUG_
    fprintf(stderr, "(2) return <false>\n");
#endif
    return false;
  }


  //  boxtest

  if( l1.p2.x > l1.p1.x ) {
    xmin = l1.p1.x;   
    xmax = l1.p2.x; 
  } else {     
    xmin = l1.p2.x; 
    xmax = l1.p1.x; 
  }

  if( l1.p2.y > l1.p1.y ) {
    ymin = l1.p1.y; 
    ymax = l1.p2.y; 
  } else {
    ymin = l1.p2.y; 
    ymax = l1.p1.y;
  }
        

  // x1 = l2.p1.x;      y1 = l2.p1.y;
  // x2 = l2.p2.x;      y2 = l2.p2.y;
    
  /*
    if( x1 > xmax  &&  x2 > xmax )  return false;
    if( x1 < xmin  &&  x2 < xmin )  return false;
    if( y1 > ymax  &&  y2 > ymax )  return false;
    if( y1 < ymin  &&  y2 < ymin )  return false;
  */


  if(l2.p1.x > xmax && l2.p2.x > xmax) {
#ifdef _DEBUG_
    fprintf(stderr, "(3) return <false>\n");
#endif
    return false;
  }

  if(l2.p1.x < xmin && l2.p2.x < xmin) {
#ifdef _DEBUG_
    fprintf(stderr, "(4) return <false>\n");
#endif
    return false;
  }

  if(l2.p1.y > ymax && l2.p2.y > ymax) {
#ifdef _DEBUG_
    fprintf(stderr, "(5) return <false>\n");
#endif
    return false;
  }

  if(l2.p1.y < ymin && l2.p2.y < ymin) {
#ifdef _DEBUG_
    fprintf(stderr, "(6) return <false>\n");
#endif
    return false;
  }

  //  get point of intersection
   
  x21 = l1.p2.x - l1.p1.x;
  x31 = l2.p1.x - l1.p1.x;
  x41 = l2.p2.x - l1.p1.x;
  x43 = l2.p2.x - l2.p1.x;

  y21 = l1.p2.y - l1.p1.y;
  y31 = l2.p1.y - l1.p1.y;
  y41 = l2.p2.y - l1.p1.y;
  y43 = l2.p2.y - l2.p1.y;
    
  d0  = x43 * y21  -  x21 * y43;
  d1  = x43 * y31  -  x31 * y43;
  d2  = x21 * y31  -  x31 * y21;


  //  check if edges are parallel

  if( fabs(d0) < CEPS && fabs(d1) > CEPS) {
#ifdef _DEBUG_
    fprintf(stderr, "(7) return <false>\n");
#endif
    return false;
  } else { 
    if( fabs(d0)< CEPS && fabs(d1) < CEPS) {
      d0 = x21 * x21  +  y21 * y21;
      d1 = x21 * x31  +  y21 * y31;
      d2 = x21 * x41  +  y21 * y41;
      s  = d1 / d0;
      t  = d2 / d0;
    } else {
      s = d1 / d0;
      t = d2 / d0;
    }
  }
        
  if( s < -CEPS  ||  s > 1.0 + CEPS ) {
#ifdef _DEBUG_
    fprintf(stderr, "(8) return <false>\n");
#endif
    return false;
  }

  if( t < -CEPS  ||  t > 1.0 + CEPS ) {
#ifdef _DEBUG_
    fprintf(stderr, "(9) return <false>\n");
#endif
    return false;
  }

#ifdef _DEBUG_
  fprintf(stderr, "<geometry::intersects()> return <true>\n");
#endif

  return true;

}

---