geometry3D/ld/line.hpp

View this file on GitHub
Last update: 2023-05-17 22:41:14+09:00
Include: #include "geometry3D/ld/line.hpp"

Depends on

Code

#pragma once

#include "base_ld.hpp"

namespace lib {

struct line {
    vec a, b;
};

vec proj(const line &l, const vec &p) {
    vec ab = l.b - l.a;
    return l.a + ab * (dot(ab, p - l.a) / norm(ab));
}

vec refl(const line &l, const vec &p) { return proj(l, p) * ld(2) - p; }

int intersection(const line &p, const line &q){
    // cross_point = alpha * p.a + (1-alpha) * p.b = beta * q.b + (1-beta) * q.a
    // alpha * vp + beta * vq = vr
    vec vp = p.a - p.b, vq = q.a - q.b, vr = q.a - p.b;
    vec cpq = cross(vp,vq), cpr = cross(vp,vr), crq = cross(vr,vq);
    // yz projection
    if (sgn(cpq.x) != 0){
        ld alpha = crq.x / cpq.x, beta = cpr.x / cpq.x;
        // cross
        if (sgn(abs(alpha * vp.x + beta * vq.x - vr.x)) == 0){
            return 0;
        }
        // nejire
        return 1;
    }
    // zx projection
    if (sgn(cpq.y) != 0){
        ld alpha = crq.y / cpq.y, beta = cpr.y / cpq.y;
        // cross
        if (sgn(abs(alpha * vp.y + beta * vq.y - vr.y)) == 0){
            return 0;
        }
        // nejire
        return 1;
    }
    // xy projection
    if (sgn(cpq.z) != 0){
        ld alpha = crq.z / cpq.z, beta = cpr.z / cpq.z;
        // cross
        if (sgn(abs(alpha * vp.z + beta * vq.z - vr.z)) == 0){
            return 0;
        }
        // nejire
        return 1;
    }
    // cpq == 0  ->  parallel
    // same
    if (sgn(abs(cross(p.a - q.a, p.b - q.a))) == 0){
        return 3;
    }
    // not same
    return 2;
}

ld dist(const line &l, const vec &p){
    return abs(p - proj(l,p));
}

vec cross_point(const line &p, const line &q){
    assert(intersection(p,q) == 0);
    // cross_point = alpha * p.a + (1-alpha) * p.b = beta * q.b + (1-beta) * q.a
    // alpha * vp + beta * vq = vr
    vec vp = p.a - p.b, vq = q.a - q.b, vr = q.a - p.b;
    vec cpq = cross(vp,vq), cpr = cross(vp,vr), crq = cross(vr,vq);
    auto res = [&](ld alpha){
        return alpha * p.a + (1-alpha) * p.b;
    };
    // yz projection
    if (sgn(cpq.x) != 0){
        ld alpha = crq.x / cpq.x, beta = cpr.x / cpq.x;
        // cross
        if (sgn(abs(alpha * vp.x + beta * vq.x - vr.x)) == 0){
            return res(alpha);
        }
    }
    // zx projection
    if (sgn(cpq.y) != 0){
        ld alpha = crq.y / cpq.y, beta = cpr.y / cpq.y;
        // cross
        if (sgn(abs(alpha * vp.y + beta * vq.y - vr.y)) == 0){
            return res(alpha);
        }
    }
    // xy projection
    if (sgn(cpq.z) != 0){
        ld alpha = crq.z / cpq.z, beta = cpr.z / cpq.z;
        // cross
        if (sgn(abs(alpha * vp.z + beta * vq.z - vr.z)) == 0){
            return res(alpha);
        }
    }
    // NOT expected
    return vec();
}


}  // namespace lib

#line 2 "geometry3D/ld/line.hpp"

#line 2 "geometry3D/ld/base_ld.hpp"

#line 2 "template/template.hpp"

#include <bits/stdc++.h>

#define rep(i, s, n) for (int i = (int)(s); i < (int)(n); i++)
#define rrep(i, s, n) for (int i = (int)(n)-1; i >= (int)(s); i--)
#define all(v) v.begin(), v.end()

using ll = long long;
using ld = long double;
using ull = unsigned long long;

template <typename T> bool chmin(T &a, const T &b) {
    if (a <= b) return false;
    a = b;
    return true;
}
template <typename T> bool chmax(T &a, const T &b) {
    if (a >= b) return false;
    a = b;
    return true;
}

namespace lib {

using namespace std;

}  // namespace lib

// using namespace lib;
#line 2 "geometry3D/base_arbitary.hpp"

#line 4 "geometry3D/base_arbitary.hpp"

namespace lib {

template<typename T>
struct Vec {
    T x, y, z;
    Vec (T _x = 0, T _y = 0, T _z = 0) : x(_x), y(_y), z(_z) {}
    Vec& operator*=(const T& a){
        x *= a;
        y *= a;
        z *= a;
        return *this;
    }
    Vec& operator/=(const T& a){
        x /= a;
        y /= a;
        z /= a;
        return *this;
    }
    Vec& operator+=(const Vec& rhs) {
        x += rhs.x;
        y += rhs.y;
        z += rhs.z;
        return *this;
    }
    Vec& operator-=(const Vec& rhs) {
        x -= rhs.x;
        y -= rhs.y;
        z -= rhs.z;
        return *this;
    }
    friend Vec operator+(const Vec& lhs, const Vec& rhs) {
        return Vec(lhs) += rhs;
    }
    friend Vec operator-(const Vec& lhs, const Vec& rhs) {
        return Vec(lhs) -= rhs;
    }
    friend Vec operator*(const Vec& lhs, const T& rhs) {
        return Vec(lhs) *= rhs;
    }
    friend Vec operator*(const T& rhs, const Vec& lhs) {
        return Vec(lhs) *= rhs;
    }
    friend Vec operator/(const Vec& lhs, const T& rhs) {
        return Vec(lhs) /= rhs;
    }
    friend std::ostream &operator<<(std::ostream &os,const Vec&r) {
        return os << "(" << r.x << "," << r.y << "," << r.z << ")"; 
    }
};

};
#line 5 "geometry3D/ld/base_ld.hpp"

namespace lib {

using vec = Vec<ld>;
const ld eps = 1e-7;

void ldout(int len = 20) { cout << fixed << setprecision(len); }

int sgn(ld a) { return (a < -eps) ? -1 : (a > eps) ? 1 : 0; }

ld dot(const vec &a, const vec &b){
    return a.x*b.x + a.y*b.y + a.z*b.z;
}

vec cross(const vec &a, const vec &b){
    return Vec(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
}

ld norm(const vec &a){
    return a.x*a.x+a.y*a.y+a.z*a.z;
}

ld abs(const vec &a){
    return sqrtl(norm(a));
}


}  // namespace lib
#line 4 "geometry3D/ld/line.hpp"

namespace lib {

struct line {
    vec a, b;
};

vec proj(const line &l, const vec &p) {
    vec ab = l.b - l.a;
    return l.a + ab * (dot(ab, p - l.a) / norm(ab));
}

vec refl(const line &l, const vec &p) { return proj(l, p) * ld(2) - p; }

int intersection(const line &p, const line &q){
    // cross_point = alpha * p.a + (1-alpha) * p.b = beta * q.b + (1-beta) * q.a
    // alpha * vp + beta * vq = vr
    vec vp = p.a - p.b, vq = q.a - q.b, vr = q.a - p.b;
    vec cpq = cross(vp,vq), cpr = cross(vp,vr), crq = cross(vr,vq);
    // yz projection
    if (sgn(cpq.x) != 0){
        ld alpha = crq.x / cpq.x, beta = cpr.x / cpq.x;
        // cross
        if (sgn(abs(alpha * vp.x + beta * vq.x - vr.x)) == 0){
            return 0;
        }
        // nejire
        return 1;
    }
    // zx projection
    if (sgn(cpq.y) != 0){
        ld alpha = crq.y / cpq.y, beta = cpr.y / cpq.y;
        // cross
        if (sgn(abs(alpha * vp.y + beta * vq.y - vr.y)) == 0){
            return 0;
        }
        // nejire
        return 1;
    }
    // xy projection
    if (sgn(cpq.z) != 0){
        ld alpha = crq.z / cpq.z, beta = cpr.z / cpq.z;
        // cross
        if (sgn(abs(alpha * vp.z + beta * vq.z - vr.z)) == 0){
            return 0;
        }
        // nejire
        return 1;
    }
    // cpq == 0  ->  parallel
    // same
    if (sgn(abs(cross(p.a - q.a, p.b - q.a))) == 0){
        return 3;
    }
    // not same
    return 2;
}

ld dist(const line &l, const vec &p){
    return abs(p - proj(l,p));
}

vec cross_point(const line &p, const line &q){
    assert(intersection(p,q) == 0);
    // cross_point = alpha * p.a + (1-alpha) * p.b = beta * q.b + (1-beta) * q.a
    // alpha * vp + beta * vq = vr
    vec vp = p.a - p.b, vq = q.a - q.b, vr = q.a - p.b;
    vec cpq = cross(vp,vq), cpr = cross(vp,vr), crq = cross(vr,vq);
    auto res = [&](ld alpha){
        return alpha * p.a + (1-alpha) * p.b;
    };
    // yz projection
    if (sgn(cpq.x) != 0){
        ld alpha = crq.x / cpq.x, beta = cpr.x / cpq.x;
        // cross
        if (sgn(abs(alpha * vp.x + beta * vq.x - vr.x)) == 0){
            return res(alpha);
        }
    }
    // zx projection
    if (sgn(cpq.y) != 0){
        ld alpha = crq.y / cpq.y, beta = cpr.y / cpq.y;
        // cross
        if (sgn(abs(alpha * vp.y + beta * vq.y - vr.y)) == 0){
            return res(alpha);
        }
    }
    // xy projection
    if (sgn(cpq.z) != 0){
        ld alpha = crq.z / cpq.z, beta = cpr.z / cpq.z;
        // cross
        if (sgn(abs(alpha * vp.z + beta * vq.z - vr.z)) == 0){
            return res(alpha);
        }
    }
    // NOT expected
    return vec();
}


}  // namespace lib