g_line.h

Go to the documentation of this file.

/* $Id $ */

#ifndef EDGE_LINE_H
#define EDGE_LINE_H

#include "m_cpp/m_vector.h" 
#include "m_cpp/m_vector_d.h" 
#include "l_cpp/l_exception.h"
#include <limits>

namespace kjb {
class Line
{
    public:

        Line() : line_params(3, 1.0) {}

        Line(const kjb::Vector & point_1, const kjb::Vector & point_2);
        
        Line(double a, double b, double c) : line_params(3, 1.0)
        {
            line_params(0) = a;
            line_params(1) = b;
            line_params(2) = c;
        }

        Line(const kjb::Vector & iparams)
        {
            if(iparams.size() != 3)
            {
                throw kjb::Illegal_argument("Line constructor, input parameter vector must have size 3");
            }
            line_params = iparams;
        }

        Line(const Line & l) : line_params(l.line_params) { }

        Line & operator=(const Line & l)
        {
            line_params = l.line_params;
            return (*this);
        }

        inline double get_a() const
        {
            return line_params(0);
        }

        inline double get_b() const
        {
            return line_params(1);
        }

        inline double get_c() const
        {
            return line_params(2);
        }
        
        inline const kjb::Vector & get_params() const
        {
            return line_params;
        }

        double compute_y_coordinate(double ix)
        {
            if( fabs(line_params(1)) < FLT_EPSILON )
            {
                throw kjb::KJB_error("Cannot compute y coordinate for a vertical line");
            }
            return (-line_params(0)*ix - line_params(2))/line_params(1);
        }

        inline void set_a(double ia)
        {
            line_params(0) = ia;
        }

        inline void set_b(double ib)
        {
            line_params(1) = ib;
        }
        
        inline void set_c(double ic)
        {
            line_params(2) = ic;
        }
        
        inline void set_line_params(const kjb::Vector & iparams)
        {
            if(iparams.size() != 3)
            {
                throw kjb::Illegal_argument("Line constructor, input parameter vector must have size 3");
            }
            line_params = iparams;
        }

        inline bool point_is_on_line(const Vector & point) const
        {
            if(point.size() < 2)
            {
                throw kjb::Illegal_argument("point_is_on_line(), point's size < 2");
            }

            return fabs(get_a()*point(0)+get_b()*point(1)+get_c()) < FLT_EPSILON;
        }

        static bool find_line_intersection(const Line & l1, const Line & l2, kjb::Vector & ints);

        double find_distance_to_point(const kjb::Vector & point) const;

        Vector project_point_onto_line(const Vector& point) const;
   
    private:

        kjb::Vector line_params;

};


bool intersect_3D_line_with_plane
(
    kjb::Vector & intersection,
    double &t,
    const kjb::Vector & point,
    const kjb::Vector & direction,
    const kjb::Vector & plane
);


template <class VectorType>
double intersect_line_with_plane_dispatch_
(
    const VectorType& line_point,
    const VectorType& line_direction,
    const VectorType& plane_point,
    const VectorType& plane_normal
)
{
    double numerator = kjb::dot((plane_point - line_point), plane_normal);
    double denominator = kjb::dot(line_direction, plane_normal);

    if(fabs(denominator) < FLT_EPSILON)
    {
        if(fabs(numerator) > FLT_EPSILON)
            return std::numeric_limits<double>::infinity();
        else
            return std::numeric_limits<double>::quiet_NaN();
    }
    return numerator/denominator;
}


inline double intersect_line_with_plane
(
    const kjb::Vector& line_point,
    const kjb::Vector& line_direction,
    const kjb::Vector& plane_point,
    const kjb::Vector& plane_normal
)
{
    return intersect_line_with_plane_dispatch_(
            line_point,
            line_direction,
            plane_point,
            plane_normal);
}

template <std::size_t D>
inline double intersect_line_with_plane
(
    const kjb::Vector_d<D>& line_point,
    const kjb::Vector_d<D>& line_direction,
    const kjb::Vector_d<D>& plane_point,
    const kjb::Vector_d<D>& plane_normal
)
{
    return intersect_line_with_plane_dispatch_(
            line_point,
            line_direction,
            plane_point,
            plane_normal);
}

inline
Vector project_point_onto_line(const Vector& A, const Vector& B, const Vector& P)
{
    IFT(A.size() == B.size() && A.size() == P.size(), Illegal_argument,
        "Cannot project point onto line: dimensions are wrong.");

    Vector AP = P - A;
    Vector AB = B - A;

    double s = dot(AP, AB) / dot(AB, AB);
    return A + s*AB;
}

} // namespace kjb

#endif