a_plane.cxx

Go to the documentation of this file.

/*
         Copyright 2010-2017 Pierre SMARS (smars@yuntech.edu.tw)
         This program is free software: you can redistribute it and/or modify
         it under the terms of the GNU General Public License as published by
         the Free Software Foundation, either version 2 of the License, or
         (at your option) any later version.
 
         This program is distributed in the hope that it will be useful,
         but WITHOUT ANY WARRANTY; without even the implied warranty of
         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
         GNU General Public License for more details.
 
         You should have received a copy of the GNU General Public License
         along with this program.  If not, see <http://www.gnu.org/licenses/>.
         */
#include "a_plane.h"
 
#include "a_debug.h"
#include <sstream>
 
//---------------------------------------------------------------------------
const std::string a_plane::help()
{
        std::ostringstream o;
        o << "********" << std::endl;
        o << "a_plane:" << std::endl;
        o << "********" << std::endl;
        o << "This is a 'plane' class" << std::endl;
        o << std::endl;
        o << "Create a plane:" << std::endl;
        o << " a_plane p p1 p2 p3" << std::endl;
        o << " (p: name of the plane)" << std::endl;
        o << " (p1, p2, p3: 3 points defining the plane)" << std::endl;
        o << std::endl;
        o << "Commands:" << std::endl;
        o << "n: returns normal to the plane." << std::endl;
        o << "dist _pt: returns distance to a point." << std::endl;
        o << "closest _pt: point on the plane closest to p." << std::endl;
        o << "intersect _seg: intersection of plane with a segment" << std::endl;
        o << "intersect _plane: intersection of plane with another plane" << std::endl;
        o << "contains _pt: check whether point is in plane" << std::endl;
        o << a_geom_base::help();
        return o.str();
}
//---------------------------------------------------------------------------
a_plane::a_plane(const a_point& p1, const a_point& p2, const a_point& p3)
{
        a_point n = cross(p2-p1,p3-p1).normalise();
        double d = n*p1;
        this->X(n.x()); 
        this->Y(n.y());
        this->Z(n.z());
        this->T(d);
}
//---------------------------------------------------------------------------
a_point a_plane::n() const
{
        return a_point(this->X(),this->Y(),this->Z());
}
//---------------------------------------------------------------------------
a_point a_plane::p0() const
{
        return this->d0()*this->n();
}
//---------------------------------------------------------------------------
double a_plane::dist(const a_point& p) const
{
        return fabs(this->T()-this->n()*p);
}
//---------------------------------------------------------------------------
a_point a_plane::closest(const a_point& p) const
{
        return p + (this->dist(p))*(this->n());
}
//---------------------------------------------------------------------------
a_point a_plane::intersect(const a_segment& s) const
{
        a_point p1 = s.p1();
        a_point n = this->n();
        a_point v = s.dir();
        double cos_angle = n*v;
        if (fabs(cos_angle) < verysmall_)
                throw a_plane::segment_parallel_error();
        double dist = this->T()-n*p1;
        v *= dist/cos_angle;
        return (p1+v);
}
//---------------------------------------------------------------------------
a_segment a_plane::intersect(const a_plane& p) const
{
        a_point n1 = this->n();
        a_point n2 = p.n();
        //check whether planes are parallel
        double cos_angle = n1*n2; 
        if (fabs(cos_angle) < verysmall_)
                throw a_plane::plane_parallel_error();
        //direction of interesction
        a_point d = cross(n1,n2);
        //direction in planes perp to direction of intersection
        a_point s1 = cross(n1,d);
        a_point s2 = cross(n2,d);
        //plane points closest to origin
        a_point p1 = this->p0();
        a_point p2 = p.p0();
        a_segment seg1(p1,p1+s1);
        a_segment seg2(p2,p2+s2);
        double k1,k2;
        a_segment seg = seg1.intersect(seg2,k1,k2);
        a_point c = seg.c();
        return a_segment(c,c+d);
}
//---------------------------------------------------------------------------
bool a_plane::contains(const a_point& p) const
{
        if (this->dist(p)<verysmall_)
                return true;
        return false;
}
//---------------------------------------------------------------------------
void a_plane::move(a_coord p)
{
        //not very clean, a_point should be a a_3dh!
        //wrong sense of rotation with a_coord, need to check other uses
        a_point t(p(0,3),p(1,3),p(2,3));
        a_point n = this->n();
        double disp = t*n;
        this->T(this->T()-disp);
        p(0,3) = p(1,3) = p(2,3) = 0.;
        a_3dh n2(n.x(),n.y(),n.z());
        a_3dh nn = p*n2;
        this->X(nn.x());
        this->Y(nn.y());
        this->Z(nn.z());
}
//***************************************************************************
//---------------------------------------------------------------------------
void a_plane::read(std::istream &in)
{
        double x,y,z,t;
        in >> x >> y >> z >> t;
        this->X(x);
        this->Y(y);
        this->Z(z);
        this->T(t);
}
//---------------------------------------------------------------------------
void a_plane::write(std::ostream &o) const
{
        o << this->X() << " ";
        o << this->Y() << " ";
        o << this->Z() << " ";
        o << this->T() << " ";
}
//-operator>>----------------------------------------------------------------
std::istream& operator>> (std::istream& i, a_plane & m)
{
        m.read(i);
        return i;
}
//-operator<<----------------------------------------------------------------
std::ostream& operator<< (std::ostream& o, const a_plane & m)
{
        m.write(o);
        return o;
}