halfspaceD2.h

Go to the documentation of this file.

#ifndef HALFSPACED2_H
#define HALFSPACED2_H

#include <cassert>
#include <sstream>
using namespace std;

#include <typedefs.h>
#include <partitionspace.h>
#include <mathlib.h>
#include <zero.h>

template< typename PT, typename PD >
class halfspaceD2 : public partitionspace<PT>
{
public:
  
  PT p0;  
  PT p1;  
  PT normal; 

  halfspaceD2() {}

  halfspaceD2( PT const & p0_, PT const & p1_ )
    { set(p0_,p1_); }
  void set( PT const & p0_, PT const & p1_ );

// Problems with template, tried explicit copy constructor.
//  halfspaceD2( halfspaceD2<PT,PD> const & h)
//    : p0(h.p0), p1(h.p1), normal(h.normal) {}

  void normalcalculate()
  {    
    PT t(p1-p0);
    normal.x = -t.y;
    normal.y = t.x;
  }

  boolc isInside( PT const & x ) const 
    { return 0 < (x.x-p0.x)*normal.x + (x.y-p0.y)*normal.y; }

  boolc isInsideOrOnBoundary( PT const & x ) const
    { return 0 < zero<PD>::val+(x.x-p0.x)*normal.x + (x.y-p0.y)*normal.y; }
    //{ return 0 < zero+(x.x-p0.x)*normal.x + (x.y-p0.y)*normal.y; }
 
  template< typename U >
  void pointOnLine
  (
    PT & x, 
    U const t
  ) const
    { x = p0 + (p1-p0)*t; }

  boolc intersectionIn_t
  ( 
    PT & t, 
    PT const & a,
    PT const & b
  ) const;

  boolc intersection
  ( 
    PT & p, 
    PT const & a,
    PT const & b
  ) const;

  boolc isOnBoundary(PT const & w) const
  { 
    PD x = (p1.y-p0.y)*(w.x-p0.x) - (p1.x-p0.x)*(w.y-p0.y);
    if ( 0 < (x + zero<PD>::val) )
    {
      if (x < zero<PD>::val)
        return true;
    }

    return false;  
  }

  boolc clip( PD & t0, PD & t1, PT const & a, PT const & m ) const;

  boolc clipNeg( PD & t0, PD & t1, PT const & a, PT const & m ) const;

  void minimizepointtoline(PD & t, PT const & w) const
  {
    PT z(p1-p0);
    t=((w.x-p0.x)*z.x+(w.y-p0.y)*z.y)/(z.x*z.x+z.y*z.y);
  }
    
  PD const distancefromhalfspace(PT const & w) const;

  operator stringc () const
    { stringstream ss; ss << p0 << " " << p1; return ss.str(); }

  stringc print() const
  { 
    stringstream ss; 
    ss << SHOW(p0) << " " << SHOW(p1) << " " << SHOW(normal); 
    return ss.str();
  }
  
};

//---------------------------------------------------------
//  Implementation.

template< typename PT, typename PD >
PD const halfspaceD2<PT,PD>::distancefromhalfspace(PT const & w) const
{ 
  PD t;  
  minimizepointtoline(t,w); 
  PT w2; 
  pointOnLine(w2,t); 
  PT w3(w2.x-w.x,w2.y-w.y);
  return w3.x*w3.x+w3.y*w3.y;
}



template< typename PT, typename PD >
void halfspaceD2<PT,PD>::set( PT const & p0_, PT const & p1_ )
{
  p0 = p0_;
  p1 = p1_;
  normalcalculate();
}

template< typename PT, typename PD >
boolc halfspaceD2<PT,PD>::intersection
( 
  PT & p, 
  PT const & a,
  PT const & b
) const
{
  PT const u(b-a);
  PT const v(p0-p1);
  PT const w(p0-a);
  // t0*u+t1*v=w
  PT t;
  bool res = d2matsolve(t,u,v,w);
  if (res==false)
    return false;

  pointOnLine(p,t.y);

  return true;
}

template< typename PT, typename PD >
boolc halfspaceD2<PT,PD>::intersectionIn_t
( 
  PT & t, 
  PT const & a,
  PT const & b
) const
{
  PT const u(b-a);
  PT const v(p0-p1);
  PT const w(p0-a);
  // t0*u+t1*v=w

  return d2matsolve(t,u,v,w);
}

template< typename PT, typename PD >
boolc halfspaceD2<PT,PD>::clipNeg
( 
  PD & t0, 
  PD & t1, 
  PT const & a, 
  PT const & m 
) const
{

  int k=0;
  if (!isInside(a+m*t0))
    k += 1;
  if (!isInside(a+m*t1))
    k += 2;

  switch (k)
  {
    case 0: return false;
    case 1: 
    {
      point2<PD> t;
      if ( solver<PD>::d2linearequ(t,m,p0-p1,p0-a) == false )
        return !isInside(a);
      t1=t[0];
      break;
    }   
    case 2:
    {
      point2<PD> t;
      if ( solver<PD>::d2linearequ(t,m,p0-p1,p0-a) == false )
        return !isInside(a);
      t0=t[0];
      break;
    }
    case 3: return true;
  }

  return true;
}


template< typename PT, typename PD >
boolc halfspaceD2<PT,PD>::clip
( 
  PD & t0, 
  PD & t1, 
  PT const & a, 
  PT const & m 
) const
{
  int k=0;
  if (isInside(a+m*t0))
    k += 1;
  if (isInside(a+m*t1))
    k += 2;

  switch (k)
  {
    case 0: return false;
    case 1: 
    {
      point2<PD> t;
      if ( solver<PD>::d2linearequ(t,m,p0-p1,p0-a) == false )
        return isInside(a);
      t1=t[0];
      break;
    }   
    case 2:
    {
      point2<PD> t;
      if ( solver<PD>::d2linearequ(t,m,p0-p1,p0-a) == false )
        return isInside(a);
      t0=t[0];
      break;
    }
    case 3: return true;
  }

  return true;
}



#endif