rLine.cpp

#include <math.h>
#include <ostream>
#include <iostream>
#include "rPoint2D.h"

#include "rPoint.h"
#include "rVector.h"
#include "rLine.h"
using std::ostream;

/* needed for MAX post.
#include <string.h>
#include <Video.h>
#include <OSUtils.h>
#include <QuickTimeComponents.h>
#include "maxincludes.h"
*/
/*
 * Class declaration^
 */
ostream& operator<<(ostream& os, rLine& p) {
  os <<"[";
  os << "<q:" <<p.q.x<<" "<<p.q.y<<" "<<p.q.z<< ">";
  os << "<n:" <<p.n.x<<" "<<p.n.y<<" "<<p.n.z<< ">";
  os <<"]";
  return os;
} /* used with cout */


rLine::rLine(rVector v)
{
  	set (v);
}
rLine::rLine(rPoint p,rVector v)
{
 	set (p,v);
}

rLine::rLine(rVector p,rVector v)
{
  set (p,v);
}

rLine::rLine(rLine &p)
{
	set (p);
}
rLine::rLine(double px,double py,double pz,double nx,double ny,double nz)
{
	set (px,py,pz,nx,ny,nz);
}
rLine& rLine::operator=(rLine p) {
	 n = p.n;
	 q = p.q; 
	 return (*this);
}
rLine::operator rVector ()
{
  return n;
}
// return the base rPoint. 
rLine::operator rPoint ()
{
  return q;
}
void rLine::set(rVector &v)
{
  //n = v;
  n.set (v.x,v.y,v.z);
  q.set (0.0,0.0,0.0);
  n.normalize ();
}
void rLine::set(rPoint &p,rVector &v)
{
  n = v;
  q = p; 
  n.normalize ();
}
void rLine::set_with_normal(rVector &p,rVector &v)
{
  n = v;
  q = p;
}

void rLine::set(rVector &p,rVector &v)
{
  //n = v;
  n.set (v.x,v.y,v.z);
  q.x = p.x;
  q.y = p.y;
  q.z = p.z;// = (rPoint)p; 
  n.normalize ();
}

void rLine::set(rLine &p)
{
  n = p.n;
  q = p.q; 
  n.normalize ();
}
void rLine::set (double px,double py,double pz,double nx,double ny,double nz)
{
  q.set (px,py,pz);
  n.set (nx,ny,nz);
  //n.normalize ();
}
// return the direction rVector.

double rLine::sub_det (rVector v2,int which)
{
  //cerr <<"n:"<<n<<"  v:"<<v2<<endl;
  switch (which) {
    case 0: return n.y*v2.x - n.x*v2.y;
    case 1: return n.z*v2.y - n.y*v2.z;
    case 2: return n.x*v2.z - n.z*v2.x;
  }
  return 0;
}
double rLine::sub_det (rVector v1,rVector v2,int which)
{
  //cerr <<"v1:"<<v1<<"  v2:"<<v2<<endl;
  switch (which) {
    case 0: return v1.y*v2.x - v1.x*v2.y;
    case 1: return v1.z*v2.y - v1.y*v2.z;
    case 2: return v1.x*v2.z - v1.z*v2.x;
  }
  return 0;
}
double rLine::intersect_param (rVector q2,rVector n2,int which)
{
  switch (which) {
    case 0: if (sub_det(n2,0)==0) return 1;
            else return ((q2.y - q.y)*n2.x + (q.x - q2.x)*n2.y)/sub_det(n2,0);
    case 1: if (sub_det(n2,1)==0) return 1;
            else return ((q2.z - q.z)*n2.y + (q.y - q2.y)*n2.z)/sub_det(n2,1);
    case 2: if (sub_det(n2,2)==0) return 1;
            else return ((q2.x - q.x)*n2.z + (q.z - q2.z)*n2.x)/sub_det(n2,2);
  }
  return 0;
}
double rLine::intersect_param (rVector q2,rVector n2,double subdet,int which)
{
  if (subdet == 0) return 0;
  //cerr <<"interparam:"<<(q2.y - q.y)<<"  "<<n2.x<<"  "<<(q.x - q2.x)<<"  "<<n2.y<<"  subdet:"<<subdet<<"  which:"<<which<<endl;
  switch (which) {
    case 0: return ((q2.y - q.y)*n2.x + (q.x - q2.x)*n2.y)/subdet;
    case 1: return ((q2.z - q.z)*n2.y + (q.y - q2.y)*n2.z)/subdet;
    case 2: return ((q2.x - q.x)*n2.z + (q.z - q2.z)*n2.x)/subdet;
  }
  return 0;
}
double rLine::second_param (rVector q2,rVector n2,double t1,int which)
{
  /* already tested to see if either n2.x or n2.y or n2.z are 0 */
  switch (which) {
    case 0: if (n2.x==0) return 0;
            return (q.x - q2.x + t1*n.x) / n2.x;
    case 1: if (n2.y==0) return 0;
            return (q.y - q2.y + t1*n.y) / n2.y;
    case 2: if (n2.z==0) return 0;
            return (q.z - q2.z + t1*n.z) / n2.z;
  }
  return 0;
}
double rLine::rev_sub_det (rVector v2,int which)
{
  //cerr <<"n:"<<n<<"  v2:"<<v2<<endl;
  switch (which) {
    case 0: return n.x*v2.y - n.y*v2.x;
    case 1: return n.y*v2.z - n.z*v2.y;
    case 2: return n.z*v2.x - n.x*v2.z;
  }
  return 0;
}
double rLine::rev_sub_det (rVector v1,rVector v2,int which)
{
  switch (which) {
    case 0: return v1.x*v2.y - v1.y*v2.x;
    case 1: return v1.y*v2.z - v1.z*v2.y;
    case 2: return v1.z*v2.x - v1.x*v2.z;
  }
  return 0;
}
double rLine::rev_intersect_param (rVector q2,rVector n2,int which)
{

  switch (which) {
    case 0: if (rev_sub_det(n2,0)==0) return 1;
            else return ((q2.x - q.x)*n2.y + (q.y - q2.y)*n2.x)/rev_sub_det(n2,0);
    case 1: if (rev_sub_det(n2,1)==0) return 1;
            else return ((q2.z - q.z)*n2.y + (q.y - q2.y)*n2.z)/rev_sub_det(n2,1);
    case 2: if (rev_sub_det(n2,2)==0) return 1;
            else return ((q2.x - q.x)*n2.z + (q.z - q2.z)*n2.x)/rev_sub_det(n2,2);
  }
  return 0;
}
double rLine::rev_intersect_param (rVector q2,rVector n2,double revsubdet,int which)
{
  if (revsubdet == 0) return 0;
  switch (which) {
    case 0: return ((q2.x - q.x)*n2.y + (q.y - q2.y)*n2.x)/revsubdet;
    case 1: return ((q2.z - q.z)*n2.y + (q.y - q2.y)*n2.z)/revsubdet;
    case 2: return ((q2.x - q.x)*n2.z + (q.z - q2.z)*n2.x)/revsubdet;
  }
  return 0;
}
double rLine::rev_second_param (rVector q2,rVector n2,double t1,int which)
{
  /* already tested to see if either n2.x or n2.y or n2.z are 0 */
  switch (which) {
    case 0: return (q.y - q2.y + t1*n.y) / n2.y;
    case 1: return (q.z - q2.z + t1*n.z) / n2.z;
    case 2: return (q.x - q2.x + t1*n.x) / n2.x;
  }
  return 0;
}


/* return the intersection of the rLine with ln.  Result is v. */
int rLine::intersect (rLine ln,rVector &v)
{
  double t1,t2;

  if (ln.q == q) {
    v = q;
    return 1;
  }
  /* if the normal rVector is 0 then the rLine is a rPoint. */
  if (ln.n.iszero ()) { /* ln is a rPoint */
    if (n.iszero ()) { /* both are rPoints */
      //cerr <<"rLine::intersect:Both rLines are rPoints\n";
      return q == ln.q;
    }
    if (intersect (ln.q,t1)) { 
      v = ln.q; 
      //cerr <<"rLine::intersect:found\n";
      return 1; 
    }
    else {
      //cerr <<"rLine::intersect:no intersection found case 1\n";
      return 0;
    }
  }
  else if (n.iszero ()) {
    if (ln.n.iszero ()) {
      return q == ln.q;
    }
    if (ln.intersect (q,t2)) { v = q; return 1; }
    else {
      //cerr <<"rLine::intersect:no intersection found case 2\n";
      return 0;
    }
  }
  rVector v2;
  double det;
  det = sub_det (ln.n,0);
  //cerr <<"det case3:"<<det<<" ln.n.x:"<<ln.n.x<<endl;
  if (det != 0 && ln.n.x != 0) {
    t1 = intersect_param (ln.q,ln.n,det,0);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,0);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 3?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = sub_det (ln.n,1);
  //cerr <<"det case4:"<<det<<" ln.n.y:"<<ln.n.y<<endl;
  if (det != 0 && ln.n.y != 0) {
    t1 = intersect_param (ln.q,ln.n,det,1);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,1);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 4?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = sub_det (ln.n,2);
  //cerr <<"det case5:"<<det<<" ln.n.z:"<<ln.n.z<<endl;
  if (det != 0 && ln.n.z != 0) {
    t1 = intersect_param (ln.q,ln.n,det,2);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,2);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 5?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  
  det = ln.sub_det (n,0);
  //cerr <<"det case6:"<<det<<" n.x:"<<n.x<<endl;
  if (det != 0 && n.x != 0) {
    t1 = ln.intersect_param (q,n,det,0);
    v = ln.solve (t1); 
    t2 = ln.second_param (q,n,t1,0);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 6?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = ln.sub_det (n,1);
  //cerr <<"det case7:"<<det<<" n.y:"<<n.y<<endl;
  if (det != 0 && n.y != 0) {
    t1 = ln.intersect_param (q,n,det,1);
    v = ln.solve (t1); 
    t2 = ln.second_param (q,n,t1,1);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 7?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = ln.sub_det (n,2);
  //cerr <<"det case8:"<<det<<" n.z:"<<n.z<<endl;
  if (det != 0 && n.z != 0) {
    t1 = ln.intersect_param (q,n,det,2);
    v = ln.solve (t1); 
    t2 = ln.second_param (q,n,t1,2);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 8?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  /* reverse variables in next 6 cases. */ 
  det = rev_sub_det (ln.n,0);
  //cerr <<"det case9:"<<det<<" ln.n.x:"<<ln.n.x<<endl;
  if (det != 0 && ln.n.x != 0) {
    t1 = rev_intersect_param (ln.q,ln.n,det,0);
    v = solve (t1); 
    t2 = rev_second_param (ln.q,ln.n,t1,0);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 9?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = rev_sub_det (ln.n,1);
  //cerr <<"det case10:"<<det<<" ln.n.y:"<<ln.n.y<<endl;
  if (det != 0 && ln.n.y != 0) {
    t1 = rev_intersect_param (ln.q,ln.n,det,1);
    v = solve (t1); 
    t2 = rev_second_param (ln.q,ln.n,t1,1);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 10?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = rev_sub_det (ln.n,2);
  //cerr <<"det case11:"<<det<<" ln.n.z:"<<ln.n.z<<endl;
  if (det != 0 && ln.n.z != 0) {
    t1 = rev_intersect_param (ln.q,ln.n,det,2);
    v = solve (t1); 
    t2 = rev_second_param (ln.q,ln.n,t1,2);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 11?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  
  det = ln.rev_sub_det (n,0);
  //cerr <<"det case12:"<<det<<" n.x:"<<n.x<<endl;
  if (det != 0 && n.x != 0) {
    t1 = ln.rev_intersect_param (q,n,det,0);
    v = ln.solve (t1); 
    t2 = ln.rev_second_param (q,n,t1,0);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 12?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = ln.rev_sub_det (n,1);
  //cerr <<"det case13:"<<det<<" n.y:"<<n.y<<endl;
  if (det != 0 && n.y != 0) {
    t1 = ln.rev_intersect_param (q,n,det,1);
    v = ln.solve (t1); 
    t2 = ln.rev_second_param (q,n,t1,1);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 13?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  det = ln.rev_sub_det (n,2);
  //cerr <<"det case14:"<<det<<" n.z:"<<n.z<<endl;
  if (det != 0 && n.z != 0) {
    t1 = ln.rev_intersect_param (q,n,det,2);
    v = ln.solve (t1); 
    t2 = ln.rev_second_param (q,n,t1,2);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 14?"<<(v2==v)<<"\n";
    return (v2==v);
  }
  else {
    //cerr <<"rLine::intersect:no intersection found, last case.\n";
    return 0;
  }
}
int rLine::intersect (rVector ln,rVector &v)
{
  return intersect (ln,v);
}
/* return the intersection of the rLine with ln.  Result is v. */
int rLine::near_intersect (rLine ln,rVector &vo1,rVector &vo2)
{
  double t1,t2;

  rVector v,v2,vd,vd2;

  if (ln.q == q) {
    v = q;
    return 1;
  }
  /* if the normal rVector is 0 then the rLine is a rPoint. */
  if (ln.n.iszero ()) { /* ln is a rPoint */
    if (n.iszero ()) { /* both are rPoints */
      //cerr <<"rLine::intersect:Both rLines are rPoints\n";
     // post ("rLine::intersect:Both rLines are rPoints");
      return q == ln.q;
    }
    if (intersect (ln.q,t1)) { 
      v = ln.q; 
      //cerr <<"rLine::intersect:found\n";
     // post ("rLine::intersect:found");
      return 1; 
    }
    else {
      //cerr <<"rLine::intersect:no intersection found case 1\n";
     // post ("rLine::intersect:no intersection found case 1");
      return 0;
    }
  }
  else if (n.iszero ()) {
    if (ln.n.iszero ()) {
      return q == ln.q;
    }
    if (ln.intersect (q,t2)) { v = q; return 1; }
    else {
      //cerr <<"rLine::intersect:no intersection found case 2\n";
     // post ("rLine::intersect:no intersection found case 2");
      return 0;
    }
  }
  double det;

  det = sub_det (ln.n,0);
  //cerr <<"det case3:"<<det<<" ln.n.x:"<<ln.n.x<<endl;
 // post ("det case3:%lf  ln.n.x:%lf  ",det,ln.n.x);
  if (det != 0 && ln.n.x != 0) {
    t1 = intersect_param (ln.q,ln.n,det,0);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,0);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 3?"<<(v2==v)<<"  t1:"<<t1<<"  t2:"<<t2<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
   //post ("rLine::intersect:case 3? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
   // post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);
    vd = (v2-v);
    //return (v2==v?1:-1);
  }
  det = sub_det (ln.n,1);
  //cerr <<"det case4:"<<det<<" ln.n.y:"<<ln.n.y<<endl;
  //post ("det case4:%lf  ln.n.y:%lf  ",det,ln.n.y);
  if (det != 0 && ln.n.y != 0) {
    t1 = intersect_param (ln.q,ln.n,det,1);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,1);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 4?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
   // post ("rLine::intersect:case 4? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
   // post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  det = sub_det (ln.n,2);
  //cerr <<"det case5:"<<det<<" ln.n.z:"<<ln.n.z<<endl;
 // post ("det case5:%lf  ln.n.z:%lf  ",det,ln.n.z);
  if (det != 0 && ln.n.z != 0) {
    t1 = intersect_param (ln.q,ln.n,det,2);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,2);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 5?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
    //post ("rLine::intersect:case 5? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
    //post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  
  det = ln.sub_det (n,0);
  //cerr <<"det case6:"<<det<<" n.x:"<<n.x<<endl;
  ///post ("det case6:%lf  n.x:%lf  ",det,n.x);
  if (det != 0 && n.x != 0) {
    t1 = ln.intersect_param (q,n,det,0);
    v = ln.solve (t1); 
    t2 = ln.second_param (q,n,t1,0);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 6?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
   // post ("rLine::intersect:case 6? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
   // post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  det = ln.sub_det (n,1);
  //cerr <<"det case7:"<<det<<" n.y:"<<n.y<<endl;
 // post ("det case7:%lf  n.y:%lf  ",det,n.y);
  if (det != 0 && n.y != 0) {
    t1 = ln.intersect_param (q,n,det,1);
    v = ln.solve (t1); 
    t2 = ln.second_param (q,n,t1,1);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 7?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
    // post ("rLine::intersect:case 7? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
    //post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  det = ln.sub_det (n,2);
  //cerr <<"det case8:"<<det<<" n.z:"<<n.z<<endl;
  //post ("det case8:%lf  n.z:%lf  ",det,n.z);
  if (det != 0 && n.z != 0) {
    t1 = ln.intersect_param (q,n,det,2);
    v = ln.solve (t1); 
    t2 = ln.second_param (q,n,t1,2);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 8?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
    //post ("rLine::intersect:case 8? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
    //post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }

  /* reverse variables in next 6 cases. */ 

  det = rev_sub_det (ln.n,0);
  //cerr <<"det case9:"<<det<<" ln.n.x:"<<ln.n.x<<endl;
  //post ("det case9:%lf  ln.n.x:%lf  ",det,ln.n.x);
  if (det != 0 && ln.n.x != 0) {
    t1 = rev_intersect_param (ln.q,ln.n,det,0);
    v = solve (t1); 
    t2 = rev_second_param (ln.q,ln.n,t1,0);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 9?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
   // post ("rLine::intersect:case 9? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
   // post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  det = rev_sub_det (ln.n,1);
  //cerr <<"det case10:"<<det<<" ln.n.y:"<<ln.n.y<<endl;
 // post ("det case10:%lf  ln.n.y:%lf  ",det,ln.n.y);
  if (det != 0 && ln.n.y != 0) {
    t1 = rev_intersect_param (ln.q,ln.n,det,1);
    v = solve (t1); 
    t2 = rev_second_param (ln.q,ln.n,t1,1);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 10?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
   // post ("rLine::intersect:case 10? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
   // post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  det = rev_sub_det (ln.n,2);
  //cerr <<"det case11:"<<det<<" ln.n.z:"<<ln.n.z<<endl;
 // post ("det case11:%lf  ln.n.z:%lf  ",det,ln.n.z);
  if (det != 0 && ln.n.z != 0) {
    t1 = rev_intersect_param (ln.q,ln.n,det,2);
    v = solve (t1); 
    t2 = rev_second_param (ln.q,ln.n,t1,2);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 11?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
    //post ("rLine::intersect:case 11? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
    //post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  
  det = ln.rev_sub_det (n,0);
  //cerr <<"det case12:"<<det<<" n.x:"<<n.x<<endl;
 // post ("det case12:%lf  n.x:%lf  ",det,n.x);
  if (det != 0 && n.x != 0) {
    t1 = ln.rev_intersect_param (q,n,det,0);
    v = ln.solve (t1); 
    t2 = ln.rev_second_param (q,n,t1,0);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 12?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
    //post ("rLine::intersect:case 12? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
    //post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  det = ln.rev_sub_det (n,1);
  //cerr <<"det case13:"<<det<<" n.y:"<<n.y<<endl;
  //post ("det case13:%lf  n.y:%lf  \n",det,n.y);
  if (det != 0 && n.y != 0) {
    t1 = ln.rev_intersect_param (q,n,det,1);
    v = ln.solve (t1); 
    t2 = ln.rev_second_param (q,n,t1,1);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 13?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
    //post ("rLine::intersect:case 13? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
    //post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }

  det = ln.rev_sub_det (n,2);
  //cerr <<"det case14:"<<det<<" n.z:"<<n.z<<endl;
  //post ("det case14:%lf  n.z:%lf  \n",det,n.z);
  if (det != 0 && n.z != 0) {
    t1 = ln.rev_intersect_param (q,n,det,2);
    v = ln.solve (t1); 
    t2 = ln.rev_second_param (q,n,t1,2);
    v2 = solve (t2);
    //cerr <<"rLine::intersect:case 14?"<<(v2==v)<<"\n";
    //cerr <<"rLine::v1:"<<v<<"  v2:"<<v2<<endl;
    //post ("rLine::intersect:case 14? %d  t1:%lf  t2:%lf ",(v2==v),t1,t2);
    //post ("rLine::v1:%lf %lf %lf   v2:%lf %lf %lf",v.x,v.y,v.z,v2.x,v2.y,v2.z);

    vd2 = (v2-v);
    //return (v2==v?1:-1);
    
  }
  if (vd2.len () < vd.len ()) { vd = vd2; vo1 = v; vo2 = v2; }
  else  { vd2 = vd; vo1 = v; vo2 = v2; }
/*
  else {
    //cerr <<"rLine::intersect:no intersection found, last case.\n";
    return 0;
  }
*/
   return (v2==v?1:-1);
}
int rLine::near_intersect (rVector ln,rVector &v,rVector &v2)
{
  return near_intersect (ln,v,v2);
}
/* return the intersection of the rLine with ln.  Result is v. */
int rLine::projected_intersect (rLine ln,rVector &v)
{
  double t1,t2;

  if (ln.q == q) {
    v = q;
    return 1;
  }
  /* if the normal rVector is 0 then the rLine is a rPoint. */
  if (ln.n.iszero ()) { /* ln is a rPoint */
    if (n.iszero ()) { /* both are rPoints */
      //cerr <<"rLine::intersect:Both rLines are rPoints\n";
      return q == ln.q;
    }
    if (projected_intersect (ln.q,t1)) { 
      v = ln.q; 
      //cerr <<"rLine::intersect:found\n";
      return 1; 
    }
    else {
      //cerr <<"rLine::intersect:no intersection found case 1\n";
      return 0;
    }
  }
  else if (n.iszero ()) {
    if (ln.n.iszero ()) {
      return q == ln.q;
    }
    if (ln.projected_intersect (q,t2)) { v = q; return 1; }
    else {
      //cerr <<"rLine::intersect:no intersection found case 2\n";
      return 0;
    }
  }
  rVector v2;
  double det;
  
  det = sub_det (ln.n,0);
  if (det != 0 && ln.n.x != 0) {
    t1 = intersect_param (ln.q,ln.n,det,0);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,0);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 3?"<<(v2==v)<<"\n";
    //post ("rLine::intersect:case 3?%d\n",(v2==v));
    return (v2==v);
  }
  det = sub_det (ln.n,1);
  if (det != 0 && ln.n.y != 0) {
    t1 = intersect_param (ln.q,ln.n,det,1);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,1);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 4?"<<(v2==v)<<"\n";
    //post ("rLine::intersect:case 4?%d\n",(v2==v));
    return (v2==v);
  }
  det = sub_det (ln.n,2);
  if (det != 0 && ln.n.z != 0) {
    t1 = intersect_param (ln.q,ln.n,det,2);
    v = solve (t1); 
    t2 = second_param (ln.q,ln.n,t1,2);
    v2 = ln.solve (t2);
    //cerr <<"rLine::intersect:case 5?"<<(v2==v)<<"\n";
    //post ("rLine::intersect:case 5?%d\n",(v2==v));
    return (v2==v);
  }
  else {
    //cerr <<"rLine::intersect:no intersection found, last case.\n";
   // post ("rLine::intersect:no intersection found, last case.\n");
    return 0;
  }
}
int rLine::projected_intersect (rVector ln,rVector &v)
{
  return projected_intersect (ln,v);
}
/*
	Scale
 */
rLine rLine::operator*(double b)
{
  rLine c;

  c.scale (b);
  return (c);
}