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

modules/transforms/vec4.asy
259 строк7 KB
Gevorkyan M. N.
init
16 июл 2026, 18:13
16 июл 2026, 18:137f048a8
struct Vec4 {
  real w, x, y, z;

  void operator init (real x, real y, real z, real w) {
    this.w = w; this.x = x; this.y = y; this.z = z;
  }

  void operator init (triple v, real w) {
    this.w = w; this.x = v.x; this.y = v.y; this.z = v.z;
  }

  void operator init (triple v) {
    this.w = 0; this.x = v.x; this.y = v.y; this.z = v.z;
  }

  void operator init (real w) {
    this.w = w; this.x = 0; this.y = 0; this.z = 0;
  }

  void operator init (real[] V) {
    this.w = V[3]; this.x = V[0]; this.y = V[1]; this.z = V[2];
  }

  // Векторная часть
  triple triple() {
    return (this.x, this.y, this.z);
  }
  triple xyz() {
    return (this.x, this.y, this.z);
  }
  // Скалярная часть
  real weight() {
    return this.w;
  }

  pair xy() {
    return (this.x, this.y);
  }

  // Конвертация проективного вектора в аффинную точку
  autounravel triple operator cast (Vec4 v) {
    if (abs(v.w) < realEpsilon) {
      return v.triple();
    } else {
      return (v.triple() / v.w);
    }
  }

  // Конвертация triple в vec4
  autounravel Vec4 operator cast (triple v) {
    return Vec4(v);
  }

  // Конвертация Vec4 в массив
  autounravel real[] operator cast (Vec4 v) {
    return new real[] {v.x, v.y, v.z, v.w};
  }

  // Конвертация массива в Vec4
  autounravel Vec4 operator cast (real[] V) {
    return Vec4(V[0], V[1], V[2], V[3]);
  }

  // Сложения
  autounravel Vec4 operator + (Vec4 v, triple u) {
    return Vec4(v.triple() + u, v.w);
  }

  autounravel Vec4 operator + (triple u, Vec4 v) {
    return Vec4(v.triple() + u, v.w);
  }

  autounravel Vec4 operator + (Vec4 u, Vec4 v) {
    return Vec4(u.triple()+v.triple(), u.w + v.w);
  }

  autounravel Vec4 operator + (Vec4 v) {
    return Vec4(v.triple(), v.w);
  }

  // Разности
  autounravel Vec4 operator - (Vec4 v, triple u) {
    return Vec4(v.triple() - u, v.w);
  }

  autounravel Vec4 operator - (triple u, Vec4 v) {
    return Vec4(u - v.triple(), -v.w);
  }

  autounravel Vec4 operator - (Vec4 u, Vec4 v) {
    return Vec4(u.triple()-v.triple(), u.w-v.w);
  }
  
  autounravel Vec4 operator - (Vec4 v) {
    return Vec4(-v.triple(), -v.w);
  }

  // Произведения
  autounravel Vec4 operator * (real a, Vec4 v) {
    return Vec4(a*v.triple(), a*v.w);
  }

  autounravel Vec4 operator * (Vec4 v, real a) {
    return Vec4(a*v.triple(), a*v.w);
  }

  // Умножение на матрицу
  autounravel Vec4 operator * (real[][] M, Vec4 v) {
    assert(
      (M.length == 4),
      "4x4 arrays only supported"
    );
    return M * (real[])v;
  }

  autounravel Vec4 operator * (Vec4 v, real[][] M) {
    assert(
      (M.length == 4),
      "4x4 arrays only supported"
    );
    return (real[])v * M;
  }

  // Скалярное произведение
  autounravel real dot(Vec4 u, Vec4 v) {
    return u.w*v.w + u.x*v.x + u.y*v.y + u.z*v.z;
  }
  // Скалярное произведение
  autounravel real dot(Vec4 u, triple v) {
    return u.x*v.x + u.y*v.y + u.z*v.z;
  }
  // Скалярное произведение
  autounravel real dot(triple u, Vec4 v) {
    return u.x*v.x + u.y*v.y + u.z*v.z;
  }

  // Векторное произведение (определено только для точек на бесконечности)
  autounravel Vec4 cross(Vec4 u, Vec4 v) {
    assert(
      (abs(u.w) <= realEpsilon) && (abs(v.w) <= realEpsilon),
      "Vectors cross product is only defined for infinity points!"
    );
    return Vec4(cross(u.triple(), v.triple()));
  }

  // Деление на действительное число
  autounravel Vec4 operator / (Vec4 v, real a) {
    assert(
      (abs(a) >= realEpsilon),
      "Division on zero!"
    );
    return Vec4(v.triple() / a, v.w / a);
  }

  // Проверка на равенство Vec4
  autounravel bool operator == (Vec4 u, Vec4 v) {
    return (u.w==v.w) && (u.x==v.x) && (u.y==v.y) && (u.z==v.z);
  }

  // real operator [] (int i) {
  //   return ((Vec4)this)[i];
  // }

  // Проверка на равенство Vec4 с учетом некоторой погрешности tol
  autounravel bool approx_equal (Vec4 u, Vec4 v, real tol=sqrt(realEpsilon)){
    bool test = true;
    test = test && (abs(u.w - v.w) <= tol);
    test = test && (abs(u.x - v.x) <= tol);
    test = test && (abs(u.y - v.y) <= tol);
    test = test && (abs(u.z - v.z) <= tol);
    return test;
  }

  // Квадрат нормы
  autounravel real abs2(Vec4 v) {
    return dot(v, v);
  }
  autounravel real norm2(Vec4 v) {
    return abs2(v);
  }
  // Норма
  autounravel real abs(Vec4 v) {
    return sqrt(abs2(v));
  }
  autounravel real norm(Vec4 v) {
    return abs(v);
  }
  // нормирование
  autounravel Vec4 unit(Vec4 v) {
    return v / abs(v);
  }
  // Проекция на единичный вектор
  Vec4 project(triple v) {
    return dot(this, unit(v))*unit(v);
  }
  // Перпендикулярная компонента относительно единичного вектора
  Vec4 perp(triple v) {
    return this - dot(this, unit(v))*unit(v);
  }
  Vec4 unitize() {
    this.x = this.x / abs(this);
    this.y = this.y / abs(this);
    this.z = this.z / abs(this);
    this.w = this.w / abs(this);
    return this;
  }
  // нормировка (приведение к евклидовом вектору/точке)
  autounravel Vec4 normalize(Vec4 v) {
    assert(
      (abs(v.w) > realEpsilon),
      "Normalization is not possible for infinity points!"
    );
    return v / v.w;
  }
  Vec4 normalize() {
    assert(
      (abs(this.w) > realEpsilon),
      "Normalization is not possible for infinity points!"
    );
    this.x /= this.w;
    this.y /= this.w;
    this.z /= this.w;
    this.w = 1;
    return this;
  }
  
  autounravel void write(Vec4 v) {
    write(v.triple(), suffix=none);
    write(", ", v.w);
  }
  // ----- Собственные движения точки и вектора ----
  Vec4 translate(triple t) {
    if (abs(this.w) <= sqrt(realEpsilon)) {
      return this;
    } else {
      this.x += t.x;
      this.y += t.y;
      this.z += t.z;
      return this;
    }
    assert(true, "Unreachable");
  }

  // Вращение вокруг оси, проходящей через начало координат
  // Vec4 rotate() {
  //   Нужно ли?
  // }
  // Вращение вокруг произвольной оси {a|ao}
  // Vec4 rotate(triple a, triple ao, real theta) {
  // Нужно ли?
  // }
}


// Заменяем triple базисные векторы на Vec4
Vec4 ex = Vec4(1, 0, 0, 0);
Vec4 ey = Vec4(0, 1, 0, 0);
Vec4 ez = Vec4(0, 0, 1, 0);
Vec4 eo = Vec4(0, 0, 0, 1);