modules/transforms/plane.asy - master - gmn/projective_geometry

modules/transforms/plane.asy
279 строк8 KB
Gevorkyan M. N.
init
16 июл 2026, 18:13
16 июл 2026, 18:137f048a8
// Зависит от 
// Vec4.asy
// lib.asy
// Line.asy
import three;
import graph3;

struct Plane {
  restricted triple n;
  restricted real d;
  restricted triple u;
  restricted triple v;
  restricted bool is_normalized = false;

  // Вычисление направляющих векторов плоскости
  void calc_uv() {
    triple[] nuv = orthogonal_triple(this.n);
    // Функция возвращает уже нормированные векторы u и v
    this.u = nuv[1];
    this.v = nuv[2];
  }

  void operator init(triple n, real d) {
    this.n = n;
    this.d = d;
    // Стоит допустить n = 0, так как это горизонт
    // надо подумать чему приравнять u и v в этом случае
    assert(
      length(n) > sqrt(realEpsilon),
      "length(n) == 0"
    );
    calc_uv();
  }
  
  // Плоскость заданная тремя точками
  void operator init(triple p1, triple p2, triple p3) {
    assert(
      abs(orient(p1, p2, p3, (0, 0, 0))) >= sqrt(realEpsilon),
      "Points are coplanar!" 
    );
    this.u = p2 - p1; 
    this.v = p3 - p1;
    this.n = cross(this.u, this.v);
    this.d = -dot(p1, this.n);
  }

  // Плоскость содержащая прямую {v|vo} и направление (ux, uy, uz | 0)
  void operator init(Line L, triple u) {
    assert(
      (abs(L.v - u) >= sqrt(realEpsilon)),
      "It is impossible to determine the plane unambiguously if L.v == u"
    );
    this.n = cross(L.v, u);
    this.d = -dot(u, L.vo);
    calc_uv();
  }

  // Плоскость содержащая прямую {v|vo} и точку (x, y, z | w)
  void operator init(Line L, Vec4 p) {
    this.n = cross(L.v, p.xyz()) + p.w*L.vo;
    this.d = -dot(p.xyz(), L.vo);
    calc_uv();
  }

  // Плоскость содержащая прямую {v|vo} и начало координат
  void operator init(Line L) {
    this.n = L.vo;
    this.d = 0;
  }

  Plane copy() {
    return Plane(this.n, this.d);
  }

  // Сложение с вектором <=> чистая трансляция
  autounravel Plane operator + (Plane P, triple t) {
    return Plane(P.n, P.d + dot(P.n, t));
  }

  Plane normalize_in_place() {
    if (this.is_normalized) {
      return this;
    } else {
      if (abs(this.n) >= realEpsilon) {
        this.d /= abs(this.n);
        this.n = unit(this.n);
      }
      else {
        this.d = 1;
      }
      this.is_normalized = true;
      return this;
    }
    assert(false, "Unreachable!");
  }

  Plane normalize() {
    if (this.is_normalized) {
      return Plane(this.n, this.d);
    } else {
      this.is_normalized = true;
      if (abs(this.n) >= realEpsilon) {
        return Plane(unit(this.n), this.d / abs(n));
      }
      else {
        return Plane((0, 0, 0), 1);
      }
    }
    assert(false, "Unreachable!");
  }

  // Расстояние от плоскости до точки
  real distance(triple p) {
    real d = dot(this.n, p) + this.d;
    return is_normalized ? d : d / abs(this.n);
  }
  
  // Расстояние от плоскости до начала координат
  real distance() {
    return -this.d / abs(this.n);
  }

  // Проекция точки на плоскость
  triple project(triple p=(0,0,0)) {
    return p - this.distance(p) * unit(this.n);
  }

  // Для рисования контура плоскости
  // по умолчанию считаем, что проекция начала координат
  // осуществляется на центр прямоугольника
  triple[] plane_as_rectangle(triple p0=(0,0,0), real su=1, real sv=1) {
    triple[] res;
    p0 = p0 + this.distance()*unit(this.n);
    res.push(p0 + 0.5(su*this.u + sv*this.v));
    res.push(p0 - 0.5(su*this.u - sv*this.v));
    res.push(p0 - 0.5(su*this.u + sv*this.v));
    res.push(p0 + 0.5(su*this.u - sv*this.v));
    return res;
  }

  /* --- Собственные движения плоскости ---
    Важно учитывать, что всякий раз, когда меняется вектор n,
    необходимо пересчитывать вектора u и v
  */
  //  Параллельный перенос на вектор t
  Plane translate(triple t) {
    // n не меняется
    this.d = d - dot(this.n, t);
    return this;
  }

  // Вращение вокруг начала координат
  Plane rotate(triple a, real theta) {
    // d не меняется
    this.n = rodrigues(p=this.n, a=a, theta=theta);
    calc_uv();
    return this;
  }

  // Винтовое движение вокруг оси вращения A
  // на угол theta, и расстояние переноса thetao
  Plane screw(Line A, real theta, real thetao) {
    Line a = A.normalize();
    triple t = sin(theta)*a.vo + (1 - cos(theta))*cross(a.v, a.vo) + thetao*a.v;
    this.n = rodrigues(p=this.n, a=a.v, theta=theta);
    this.d = this.d - dot(this.n, t);
    calc_uv();
    return this;
  }

  // --- Несобственные движения плоскости ---
  // Отражение относительно другой плоскости
  Plane reflect(Plane P) {
    Plane mirror = P.normalize();
    this.d = 2mirror.d*dot(this.n, mirror.n) - this.d;
    this.n = -reflect(mirror.n, this.n);
    calc_uv();
    return this;
  }

  // Матрица отражения точки относительно самой плоскости
  real[][] reflect() {
    return reflect(this.n, this.d);
  }

  // Зеркальный поворот (зеркальное вращение)
  Plane rotoreflection() {
    // TODO
    return this;
  }

  // Скользящее отражение (скользящая симметрия)
  Plane transflection() {
    // TODO
    return this;
  }

  // Проверка на равенство двух плоскостей с учетом некоторой погрешности величины tol
  autounravel bool approx_equal (Plane P1, Plane P2, real tol=sqrt(realEpsilon)){
    bool test = true;
    test = test && (abs(P1.n - P2.n) <= tol);
    test = test && (abs(P1.d - P2.d) <= tol);
    return test;
  }

  autounravel void write(Plane P) {
    write(P.n, suffix=none);
    write(", ", P.d);
  }
}

// ------------------------
// Вспомогательные функции
// ------------------------
// Базисные плоскости
Plane Oxy() {return Plane(Z, 0);}
Plane Oxz() {return Plane(Y, 0);}
Plane Oyz() {return Plane(X, 0);}

// Представление плоскости в виде вектора
Vec4 operator cast(Plane P) {
  return Vec4(P.n, P.d);
}

// Точка пересечения плоскости и прямой
Vec4 meet(Plane P, Line L) {
  return Vec4(
    cross(L.vo, P.n) + P.d * L.v,
    -dot(P.n, L.v)
  );
}

Vec4 meet(Line L, Plane P) {
  return meet(P, L);
}

// Плоскость, содержащая точку и прямую
Plane join(Line L, Vec4 p) {
  return Plane(L, p);
}

// Прямая пересечения двух плоскостей
Line meet(Plane P1, Plane P2) {
  return Line(
    v =cross(P1.n, P2.n),
    vo=P1.d*P2.n-P2.d*P1.n 
  );
}

// Точка пересечения трех плоскостей
Vec4 meet(Plane P1, Plane P2, Plane P3) {
  return Vec4(
    v=P1.d*cross(P3.n, P2.n) + P2.d*cross(P1.n, P3.n) + P3.d*cross(P2.n, P1.n),
    w=dot(cross(P1.n, P2.n), P3.n)
  );
}

// Рисуем ортогональную тройку векторов на поверхности плоскости
void draw_nuv(Plane P, real arrow_size=5, pen[] pens=new pen[] {deepred, deepgreen, deepblue}) {
  triple O = (0, 0, 0);
  triple p0 = P.project(O);
  // draw(O--p0, arrow=Arrow3(arrow_size));
  draw(p0--p0+P.u, p=pens[0], arrow=Arrow3(arrow_size));
  draw(p0--p0+P.v, p=pens[1], arrow=Arrow3(arrow_size));
  draw(p0--p0+P.n, p=pens[2], arrow=Arrow3(arrow_size));
}

// Перегрузка функции рисования поверхности для плоскости
void draw(
    picture pic=currentpicture, Plane P,
    triple p0=(0, 0, 0), real su=1, real sv=1,
    int nu=1, int nv=1, material surfacepen=currentpen,
    pen meshpen=nullpen, light light=currentlight,
    light meshlight=nolight, string name="",
    render render=defaultrender
  ) {
  surface PL = surface(graph(P.plane_as_rectangle(p0=p0, su, sv), join=Straight)--cycle);
  draw(pic=pic, s=PL, nu=nu, nv=nv, surfacepen=surfacepen, meshpen=meshpen, light=light, meshlight=meshlight, name=name, render=render);
}