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

modules/transforms/lib.asy
156 строк5 KB
Gevorkyan M. N.
init
16 июл 2026, 18:13
16 июл 2026, 18:137f048a8
// Алгоритм формирования ортогональной тройки векторов по заданному вектору
// Marschner S., Shirley P. Fundamentals of Computer Graphics, p 29
triple[] orthogonal_triple(triple u) {
    triple e, v, w;
    u = unit(u);
    if ( u.x <= u.y && u.x <= u.z) { e = (1, 0, 0); }
    else if (u.y <= u.x && u.y <= u.z) { e = (0, 1, 0); }
    else { e = (0, 0, 1); }
    v = unit(cross(u, e));
    w = unit(cross(u, v));
    return new triple[] {u, v, w};
}

triple operator cast(real[] v) {
  if (abs(v[3]) < 1e-10) {
    return (v[0], v[1], v[2]);
  }
  else {
    return (v[0], v[1], v[2])/v[3];
  }
}

// Сравнение двух векторов с погрешностью
bool approx_equal (triple u, triple v, real tol=sqrt(realEpsilon)){
    bool test = true;
    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;
}

// Кососимметричная матрица из вектора
real[][] skew(triple a) {
  return new real[][] {
    {0, -a.z, a.y},
    {a.z, 0, -a.x},
    {-a.y, a.x, 0}
  };
}

// Тензорное произведение векторов
// или операторное произведение векторов, или диада Гиббса
real[][] operator_product(triple a, triple b) {
  return new real[][] {
    {a.x*b.x, a.x*b.y, a.x*b.z},
    {a.y*b.x, a.y*b.y, a.y*b.z},
    {a.z*b.x, a.z*b.y, a.z*b.z}
  };
}

// Формулы Родрига (вращения) в векторном виде

// Основная формула Родрига
triple rodrigues(triple p, triple a, real theta) {
  a = unit(a);
  return cos(theta)*p + (1-cos(theta))*dot(a,p)*a  + sin(theta)*cross(a, p);
}

triple[] rodrigues(triple[] Ps, triple a, real theta) {
  triple[] res;
  for (triple p : Ps) {
    res.push(rodrigues(p=p, a=a, theta=theta));
  }
  return res;
}

// Формула Родрига для вектора конечного поворота
triple rodrigues(triple p, triple theta) {
  return p + cross(2theta/(1+abs2(theta)), p + cross(theta, p));
}
triple[] rodrigues(triple[] Ps, triple theta) {
  triple[] res;
  for (triple p : Ps) {
    res.push(rodrigues(p=p, theta=theta));
  }
  return res;
}

// Формула Родрига-Гамильтона
triple rodrigues_hamilton(triple p, triple lambda, real lambda0) {
  return p + 2lambda0*cross(lambda, p) + 2cross(lambda, cross(lambda, p));
}

triple[] rodrigues_hamilton(triple[] Ps, triple lambda, real lambda0) {
  triple[] res;
  for (triple p : Ps) {
    res.push(rodrigues_hamilton(p=p, lambda=lambda, lambda0=lambda0));
  }
  return res;
}

// Формула Родрига для поворота точки вокруг произвольной оси {a|ao}
triple rodrigues(triple p, triple a, triple ao, real theta) {
  ao = ao / abs(a);
  a = unit(a);
  return rodrigues(p, a, theta) + sin(theta)*ao + (1-cos(theta))*cross(a, ao);
}

triple[] rodrigues(triple[] Ps, triple a, triple ao, real theta) {
  triple[] res;
  for (triple p : Ps) {
    res.push(rodrigues(p=p, a=a, ao=ao, theta=theta));
  }
  return res;
}

// Формулы Родрига в матричном виде
real[][] rodrigues(triple keyword a, real keyword theta) {
  real[][] I = identity(3);
  a = unit(a);
  real[][] A = skew(a);
  return I + sin(theta) * A + (1-cos(theta)) * A * A;
}

// Формулы Родрига в матричном виде для поворота вокруг оси {a|ao}
real[][] rodrigues(real keyword theta, triple keyword a, triple keyword ao) {
  real[][] R = rodrigues(theta=theta, a=a);
  ao = ao / abs(a);
  a = unit(a);
  triple t = sin(theta) * ao + (1 - cos(theta))*cross(a, ao);
  return new real[][] {
    {R[0][0], R[0][1], R[0][2], t.x},
    {R[1][0], R[1][1], R[1][2], t.y},
    {R[2][0], R[2][1], R[2][2], t.z},
    {0, 0, 0, 1}
  };
}

// Винтовое преобразование точки общего вида реализованно в структуре line

// Формулы отражения

// для плоскости, проходящей через начало координат
triple reflect(triple n, triple p) {
  return p - 2dot(n, p)*n;
}

real[][] reflect(triple n) {
  return identity(3) - 2operator_product(n, n);
}

// Для произвольной плоскости
triple reflect(triple n, real d, triple p) {
  return p - 2(dot(n, p) + d)*n;
}

real[][] reflect(triple n, real d) {
  real[][] Rf = reflect(n);
  triple t = -2d*n;
  return new real[][] {
    {Rf[0][0], Rf[0][1], Rf[0][2], t.x},
    {Rf[1][0], Rf[1][1], Rf[1][2], t.y},
    {Rf[2][0], Rf[2][1], Rf[2][2], t.z},
    {0, 0, 0, 1}
  };
}