line.h

/********************************************************************
 * Copyright (C) 2015 Liangliang Nan <liangliang.nan@gmail.com>
 * https://3d.bk.tudelft.nl/liangliang/
 *
 * This file is part of Easy3D. If it is useful in your research/work,
 * I would be grateful if you show your appreciation by citing it:
 * ------------------------------------------------------------------
 *      Liangliang Nan.
 *      Easy3D: a lightweight, easy-to-use, and efficient C++ library
 *      for processing and rendering 3D data.
 *      Journal of Open Source Software, 6(64), 3255, 2021.
 * ------------------------------------------------------------------
 *
 * Easy3D is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License Version 3
 * as published by the Free Software Foundation.
 *
 * Easy3D 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/>.
 ********************************************************************/
 
#ifndef EASY3D_CORE_LINE_H
#define EASY3D_CORE_LINE_H
 
#include <easy3d/core/vec.h>
#include <easy3d/core/constant.h>
#include <easy3d/util/logging.h>
 
 
namespace easy3d {
 

    template <int DIM, typename FT>
    class GenericLine {
    public:
        typedef Vec<DIM, FT> Point;                
        typedef Vec<DIM, FT> Vector;          
        typedef GenericLine<DIM, FT> thisclass; 
 
    public:
        static GenericLine from_point_and_direction(const Point& p, const Vector& dir) { return GenericLine(p, dir); }
        static GenericLine from_two_points(const Point& p, const Point& q)   { return GenericLine(p, q - p); }

        GenericLine() = default;

        void set(const Point& p, const Vector& dir) {
            p_ = p;
            dir_ = normalize(dir);
        }

        const Vector& direction() const { return dir_; }

        const Point& point() const { return p_; }

        Point projection(const Point &p) const { return p_ + dir_ * dot(p - p_, dir_); }

        FT squared_distance(const Point &p) const { return length2(projection(p) - p); }

        bool feet(const thisclass& other, Point& p1, Point& p2) const;
 
    private:  // Ambiguities exist for this one.
        GenericLine(const Point & p, const Vector & dir);
 
    private:
        Point p_;      
        Vector dir_;  
    };
 
 
    template <int DIM, typename FT>
    GenericLine<DIM, FT>::GenericLine(const Point & p, const Vector & dir) : p_(p) {
      dir_ = normalize(dir);
      LOG_IF(length(dir_) < 1e-15, ERROR)
              << "degenerate line constructed from point (" << p << ") and direction (" << dir << ")";
    }
 
 
    template <int DIM, typename FT>
    bool GenericLine<DIM, FT>::feet(const thisclass& other, Point& p1, Point& p2) const {
        // See the derivation: https://www.jianshu.com/p/34a7c4e1f3f5
        FT a = dot(dir_, other.direction());
        FT b = dot(dir_, dir_);
        FT c = dot(other.direction(), other.direction());
        if (std::abs(a * a - b * c) < epsilon<FT>()) // the two lines are colinear or parallel
            return false;
        FT d = dot(other.point() - p_, dir_);
        FT e = dot(other.point() - p_, other.direction());
        FT t1 = 0;
        FT t2 = 0;
        if (a == 0) {
            t1 = d / b;
            t2 = -e / c;
        }
        else {
            t1 = (a * e - c * d) / (a * a - b * c);
            t2 = b / a * t1 - d / a;
        }
        p1 = p_ + t1 * dir_;
        p2 = other.point() + t2 * other.direction();
        return true;
    }
 

    template <int DIM, typename FT>
    std::ostream& operator<<(std::ostream& os, const GenericLine<DIM, FT>& line) {
        return os << line.point() << " " << line.direction();
    }

    template <int DIM, typename FT>
    std::istream& operator>>(std::istream& is, GenericLine<DIM, FT>& line) {
        typename GenericLine<DIM, FT>::Point  p;
        typename GenericLine<DIM, FT>::Vector dir;
        is >> p >> dir;
        line.set(p, dir);
        return is;
    }
 
}
 
#endif  // EASY3D_CORE_LINE_H