a_triangle.cxx

Go to the documentation of this file.

/*
         Copyright 2000-19 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_triangle.h"
#include "a_plane.h"
#include "a_segment.h"
 
#include "a_debug.h"
#include <sstream>
 
//---------------------------------------------------------------------------
const std::string a_triangle::help()
{
        std::ostringstream o;
        o << "***********" << std::endl;
        o << "a_triangle:" << std::endl;
        o << "***********" << std::endl;
        o << "This is a 'triangle' class" << std::endl;
        o << std::endl;
        o << "Create a triangle:" << std::endl;
        o << " a_triangle a" << std::endl;
        o << " a_triangle a _p1 _p2 _p3" << std::endl;
        o << std::endl;
        o << "Commands:" << std::endl;
        o << "p1, p2, p3: get or set apices" << std::endl;
        o << "s1, s2, s3: returns sides of triangle as segments" << std::endl;
        o << "s _ref: returns sides of triangle as segments" << std::endl;
        o << "c: get centre" << std::endl;
        o << "translate _x _y _z: translate." << std::endl;
        o << "rotate _xaxis _yaxis _zaxis: rotate." << std::endl;
        o << "invert: inverts orientation (p2 and p3)." << std::endl;
        o << "normal: get normal vector" << std::endl;
        o << "contains _pt: check whether contains point" << std::endl;
        o << "closest _pt: triangle point closest to p" << std::endl;
        o << "closestp _pt: point of plane defined by triangle closest to p" << std::endl;
        o << "dist _pt: distance between p and triangle" << std::endl;
        o << "distp _pt: distance between p and plane defined by triangle" << std::endl;
        o << "intersect _seg: intersection of triangle with a segment" << std::endl;
        o << "min_edge: length of shortest edge" << std::endl;
        o << "max_edge: length of longest edge" << std::endl;
        o << "quality: ratio between the shortest and longest edge" << std::endl;
        o << a_geom_base::help();
        return o.str();
}
//---------------------------------------------------------------------------
bool a_triangle::operator==(const a_triangle& p)
{
        return ((p1_ == p.p1())&&(p2_ == p.p2())&&(p3_ == p.p3()));
}
//---------------------------------------------------------------------------
a_triangle & a_triangle::operator=(const a_triangle& p)
{
        p1_ = p.p1();
        p2_ = p.p2();
        p3_ = p.p3();
        return *this;
}
//---------------------------------------------------------------------------
a_segment a_triangle::s(const int ref) const
{
        if (ref==0) return this->s1();
        else if (ref==1) return this->s2();
        return this->s3();
}
//-rotate--------------------------------------------------------------------
a_triangle& a_triangle::translate(double x, double y, double z)
{
        p1_.translate(x,y,z);
        p2_.translate(x,y,z);
        p3_.translate(x,y,z);
        return *this;
}
//-rotate--------------------------------------------------------------------
a_triangle& a_triangle::rotate(const a_point& x_axis, const a_point& y_axis, const a_point& z_axis)
{
        p1_.rotate(x_axis, y_axis, z_axis);
        p2_.rotate(x_axis, y_axis, z_axis);
        p2_.rotate(x_axis, y_axis, z_axis);
        return *this;
}
//---------------------------------------------------------------------------
a_triangle& a_triangle::invert()
{
        a_point v = p2_;
        p2_ = p3_;
        p3_ = v;
        return *this;
}
//---------------------------------------------------------------------------
a_point a_triangle::normal() const
{
        return cross(p2_-p1_,p3_-p1_).normalise();
}
//---------------------------------------------------------------------------
double a_triangle::S() const
{
        return cross(p2_-p1_,p3_-p1_).norm()/2.;
}
//---------------------------------------------------------------------------
a_triangle a_triangle::shape() const
{
        a_point c = this->c();
        return a_triangle(p1_-c,p2_-c,p3_-c);
}
//---------------------------------------------------------------------------
a_mat_sq a_triangle::coord() const
{
        a_mat_sq a(3);
        for (short i = 0; i<3; i++)
        {
                a_point p = this->p(i);
                for (short j = 0; j<3; j++)
                        a(i,j) = p[j];
        }
        return a;
}
//---------------------------------------------------------------------------
a_mat_sq a_triangle::inertia() const
{
        a_mat_sq a(3);
        a_triangle tri = this->shape();
        a_point c = this->c();
        double S = this->S();
        //matrix of inertia of unit triangle: [(-1/3,-1/3) (2/3,-1/3) (-1/3,2/3)]
        a_mat_sq s(3);
        for (short i = 0; i<2; i++)
        {
                s(i,i) = 1./36.;
                s(i,(i+1)%2) = -1./72.;
        }
        //matrix of transformation
        a_mat_sq v(3);
        a_point p1 = tri.p1();
        a_point p2 = tri.p2();
        for (short i = 0; i<2; i++)
        {
                v(i,0) = -2.*p1[i]-p2[i];
                v(i,1) = -p1[i]+p2[i];
                v(i,2) = 0.;
        }
        //its transpose
        a_mat_sq vt(3);
        vt = v.transpose();
        a_mat_sq cm(3);
        cm = v*s;
        cm *= vt;
        cm *= 2.*S;
        //translation
        a_mat_sq im(3);
        im = c.inertia();
        im *= S;
        a = im+cm;
        return a;
}
//---------------------------------------------------------------------------
bool a_triangle::contains(const a_point& p) const
{
        if (this->distp(p)>verysmall_)
                return false;
        a_point v1 = this->p3()-this->p1();
        a_point v2 = this->p2()-this->p1();
        a_point vp = p-this->p1();
        double a11 = v1*v1;
        double a12 = v1*v2;
        double a22 = v2*v2;
        double p1 = vp*v1;
        double p2 = vp*v2;
        double det = a11*a22-a12*a12;
        double u = (a22*p1-a12*p2)/det;
        double v = (-a12*p1+a11*p2)/det;
        if (u<0.)
                return false;
        if (v<0.)
                return false;
        if (u+v>1.)
                return false;
        return true;
}
//---------------------------------------------------------------------------
a_point a_triangle::closestp(const a_point& p) const
{
        a_plane a(this->p1(),this->p2(),this->p3());
        return a.closest(p);
}
//---------------------------------------------------------------------------
a_point a_triangle::closest(const a_point& p) const
{
        a_point pt = this->closestp(p);
        if (this->contains(pt))
                return pt;
        a_segment s1(this->p1(),this->p2());
        a_segment s2(this->p2(),this->p3());
        a_segment s3(this->p3(),this->p1());
        a_point p1 = s1.closest(p);
        a_point p2 = s2.closest(p);
        a_point p3 = s3.closest(p);
        double d1 = (p-p1).norm();
        double d2 = (p-p2).norm();
        double d3 = (p-p3).norm();
        if ((d1<d2)&&(d1<d3))
                return p1;
        if ((d2<d1)&&(d2<d3))
                return p2;
        return p3;
}
//---------------------------------------------------------------------------
double a_triangle::dist(const a_point& p) const
{
        return (p-this->closest(p)).norm();
}
//---------------------------------------------------------------------------
double a_triangle::distp(const a_point& p) const
{
        return (p-this->closestp(p)).norm();
}
//---------------------------------------------------------------------------
a_segment a_triangle::intersect(const a_segment& s) const
{
        a_point p1 = s.p1();
        a_point p2 = s.p2();
        double d1 = this->distp(p1);
        double d2 = this->distp(p2);
        if ((d1>verysmall_)||(d2>verysmall_))
        {
                throw a_segment::null_segment_error();
        }
        bool t1 = this->contains(p1); 
        bool t2 = this->contains(p2);
        if (t1&&t2)
                return s;
        std::vector<a_point> list;
        a_segment s1(p1_,p2_);
        a_segment s2(p2_,p3_);
        a_segment s3(p3_,p1_);
        double m1,m2;
        a_segment i1 = s1.shortest(s,m1,m2);
        a_segment i2 = s2.shortest(s,m1,m2);
        a_segment i3 = s3.shortest(s,m1,m2);
        if (i1.length()<verysmall_)
                list.push_back(i1.c());
        if (i2.length()<verysmall_)
                list.push_back(i2.c());
        if (i3.length()<verysmall_)
                list.push_back(i3.c());
        if (!t1&&!t2)
        {
                if (list.size()==2)
                        return a_segment(list[0],list[1]);
                else
                        throw a_segment::null_segment_error();//no intersection
        }
        if (t1)
                return a_segment(p1,list[0]);
        else
                return a_segment(p2,list[0]);
}
//---------------------------------------------------------------------------
double a_triangle::min_edge() const
{
        double m = this->s1().length();
        double d2 = this->s2().length();
        double d3 = this->s3().length();
        if (d2<m)
                m = d2;
        if (d3<m)
                m = d3;
        return m;
}
//---------------------------------------------------------------------------
double a_triangle::max_edge() const
{
        double M = this->s1().length();
        double d2 = this->s2().length();
        double d3 = this->s3().length();
        if (d2>M)
                M = d2;
        if (d3>M)
                M = d3;
        return M;
}
//---------------------------------------------------------------------------
double a_triangle::quality() const
{
        double m = this->s1().length();
        double d2 = this->s2().length();
        double d3 = this->s3().length();
        double M = m;
        if (d2<m)
                m = d2;
        else
                M = d2;
        if (d3<m)
                m = d3;
        if (d3>M)
                M = d3;
        if (M==0)
                return 0.;
        else
                return m/M;
}
//***************************************************************************
//---------------------------------------------------------------------------
void a_triangle::read(std::istream &in)
{
        in >> p1_;
        in >> p2_;
        in >> p3_;
}
//---------------------------------------------------------------------------
void a_triangle::write(std::ostream &o) const
{
        o << this->p1() << " ";
        o << this->p2() << " ";
        o << this->p3() << " ";
}