Как найти точку в прямоугольнике, которая не изменит своего положения при его уменьшении и повороте?
Решаю задачу из Тинькофф.Контест. Условие следующее:
Миша сидел в переговорке и решил нарисовать ее план на листке бумаги. Когда он закончил рисовать план переговорки, он положил лист с планом на пол переговорки. Теперь ему стало интересно: а есть ли на плане точка, которая лежит ровно на том месте пола, за которое она отвечает? Помогите ему найти эту точку.
Формат входных данных:
Переговорка и ее план имеют форму прямоугольника. Первая строка входного файла содержит два вещественных числа: координаты X и Y переговорки. Координаты углов переговорки равны (0,0),(X,0),(X,Y),(0,Y).
Вторая строка - восемь вещественных чисел, описывающих положение углов плана переговорки.Пример данных:
Ввод:
10 5
3.0 2.5 1.0 2.5 1.0 1.5 3.0 1.5
Вывод: 2.5000 2.0833
Я решил задачу для случаев, когда ось X плана переговорки параллельна оси Y или X самой переговорки. Т.е в этих случаях:
Вот мой код:
public static (double x, double y) Task8(double x, double y, double[] corners)
{
// заношу координаты планерки в массив, чтоб было удобнее обращаться к углам
// [0, 0] = (0, 0), [1, 0] = (x, 0), [0, 1] = (0, y), [1, 1] = (x, y)
(double x, double y)[,] planCornersCoordinates = new (double x, double y)[2, 2];
planCornersCoordinates[0, 0] = (corners[0], corners[1]);
planCornersCoordinates[1, 0] = (corners[2], corners[3]);
planCornersCoordinates[1, 1] = (corners[4], corners[5]);
planCornersCoordinates[0, 1] = (corners[6], corners[7]);
// Вычисляю масштаб
double planXLenght = Math.Abs(planCornersCoordinates[0, 0].x - planCornersCoordinates[1, 1].x);
double planYLenght = Math.Abs(planCornersCoordinates[0, 0].y - planCornersCoordinates[1, 1].y);
double planSquare = planXLenght * planYLenght;
double square = x * y;
double scale = Math.Sqrt(square / planSquare);
// Если ось X плана планерки параллельна оси X самой планерки, то просто решаю два уравнения с одной переменной
if (planCornersCoordinates[0, 0].y == planCornersCoordinates[1, 0].y)
{
var signX = -Math.Sign(planCornersCoordinates[1, 0].x - planCornersCoordinates[0, 0].x);
var signY = -Math.Sign(planCornersCoordinates[0, 1].y - planCornersCoordinates[0, 0].y);
x = (planCornersCoordinates[0, 0].x * scale) / (scale + signX);
y = (planCornersCoordinates[0, 0].y * scale) / (scale + signY);
}
// Если ось X плана планерки параллельна оси Y самой планерки, то решаю систему из двух уравнений
// с двумя переменными методом Крамера
if (planCornersCoordinates[0, 0].x == planCornersCoordinates[1, 0].x)
{
var signX = Math.Sign(planCornersCoordinates[1, 0].y - planCornersCoordinates[0, 0].y);
var signY = Math.Sign(planCornersCoordinates[0, 1].x - planCornersCoordinates[0, 0].x);
var determinant = 1 / scale * signX * signY * 1 / scale - 1;
var determinantX = -planCornersCoordinates[0, 0].y * 1 / scale * signY - planCornersCoordinates[0, 0].x;
x = determinantX / determinant;
var determinantY = -planCornersCoordinates[0, 0].x * 1 / scale * signX - planCornersCoordinates[0, 1].y;
y = determinantY / determinant;
}
return (x, y);
}
Но я не могу решить задачу, когда план повернут на какой-нибудь угол. Примерно такая ситуация:
Я пытаюсь составить уравнение с использованием косинуса этого угла, но так и не понял, как это нужно сделать.
Я создал репозиторий с полным условием задачи, кодом и юнит-тестами: https://github.com/ByMyKolaps/Task8_Tinkoff
Ответы (6 шт):
Если сжатие пропорционально и порядок углов плана соответствует порядку 00-X0-XY-0Y, то можно вычислить
scale = len(edge1) / X
rot = atan2(y1-y0, x1-x0)
shift = x0, y0
и составить матрицу аффинного преобразования M как произведение трех матриц выше (в таком порядке)
sc 0 0
mscl = 0 sc 0
0 0 1
cos(rot) -sin(rot) 0
mrot = sin(rot) cos(rot) 0
0 0 1
1 0 0
mshft = 0 1 0
x0 y0 1
Затем решить систему уравнений
(tx, ty, 1) * [M] = (tx, ty, 1)
Как решить без матриц.
Пусть в системе координат сторон прямоугольника координаты некоторой точки (a,b) - то есть по оси Ox координата точки a*AB, а по оси Oy - b*AD. Тогда у точки, в которую переходит данна, координаты внутри прямоугольника будут тоже (a,b): x' = A'_x + a*A'B', y' = A'_y + b*A'D'
Для поиска неподвижной точки достаточно приравнять x==x' && y==y'. Система уравнений относительно a и b на рисунке выше.
Прошу простить, в C# плаваю, писал последний раз лет 10 назад, поэтому я вам дам решение на Пайтоне, хорошо?
Для начала класс Матрица-два-на-два, чтобы решать по формуле Крамера
class Matrix2x2:
def __init__(self, rows) -> None:
self.rows = rows
def det(self):
"Детерминант"
return self[0,0]*self[1,1]-self[0,1]*self[1,0]
def __getitem__(self, idx):
"Функция индексации: позволяет писать `m[0,1]` для коэффициента в первой строке и втором столбце."
if isinstance(idx, tuple):
return self.rows[idx[0]][idx[1]]
raise ValueError("Not a tuple: ", idx)
def replace_col(self, col, a1, a2):
"Меняет колонку с номером col на столбец [a1, a2] (транспонированный)"
if col == 0:
return Matrix2x2([[a1, self[0,1]],[a2, self[1,1]]])
elif col == 1:
return Matrix2x2([[self[0,0], a1],[self[1,0], a2]])
else:
raise ValueError("Out of bounds: ", col)
def solve_eq(self, x, y):
"Решает систему линейных уравнений с правой частью [x, y]"
d0 = self.det()
m1 = self.replace_col(0, x,y)
m2 = self.replace_col(1, x,y)
return m1.det()/d0, m2.det()/d0
def __repr__(self):
"ToString"
return f"{self.rows}"
Теперь решатель и ваши тестовые случаи:
tests = [
[(10, 5), ( 4.5, 2.0, 4.5, 4.0, 5.5, 4.0, 5.5, 2.0 ), (5.104166, 3.020833)],
[(10, 5), ( 4.5, 4.0, 4.5, 2.0, 5.5, 2.0, 5.5, 4.0 ), (5.09615, 2.98076)],
[(10, 5), ( 5.5, 4.0, 5.5, 2.0, 4.5, 2.0, 4.5, 4.0 ), (4.89583, 3.020833)],
[(10, 5), ( 5.5, 2.0, 5.5, 4.0, 4.5, 4.0, 4.5, 2.0 ), (4.90384, 2.98076)],
[(10, 5), ( 1.0, 1.5, 3.0, 1.5, 3.0, 2.5, 1.0, 2.5 ), (1.25, 1.875)],
[(10, 5), ( 3.0, 1.5, 1.0, 1.5, 1.0, 2.5, 3.0, 2.5 ), (2.5, 1.875)],
[(10, 5), ( 3.0, 2.5, 1.0, 2.5, 1.0, 1.5, 3.0, 1.5 ), (2.5, 12.5 / 6)],
[(10, 5), ( 1.0, 2.5, 3.0, 2.5, 3.0, 1.5, 1.0, 1.5 ), (1.25, 12.5 / 6)]
]
def solve(test_case):
(A,D),(Ax,Ay,Bx,By,_,_,Dx,Dy),_ = test_case
m = Matrix2x2([
[A-Bx+Ax, -Dx+Ax],
[-By+Ay, D-Dy+Ay]
])
x,y = m.solve_eq(Ax,Ay)
return x*A, y*D
def test(test_case):
print(solve(test_case),test_case[-1])
for tc in tests:
test(tc)
Результат. В каждой строке первая пара - найденное решение, вторая пара ожидаемое значение.
(5.104166666666666, 3.020833333333333) (5.104166, 3.020833)
(5.096153846153846, 2.980769230769231) (5.09615, 2.98076)
(4.895833333333333, 3.020833333333333) (4.89583, 3.020833)
(4.903846153846153, 2.980769230769231) (4.90384, 2.98076)
(1.25, 1.875) (1.25, 1.875)
(2.5, 1.875) (2.5, 1.875)
(2.5, 2.0833333333333335) (2.5, 2.0833333333333335)
(1.25, 2.0833333333333335) (1.25, 2.0833333333333335)
Вот решение вашей задачи:
public static (double x, double y) Task8(double x, double y, double[] corners) {
double xx = 0;
double yy = 0;
for (int i = 0; i < 100; ++i) {
double alpha = xx / x;
double beta = yy / y;
xx = (1 - alpha) * ((1 - beta) * corners[0] + beta * corners[6]) +
alpha * ((1 - beta) * corners[2] + beta * corners[4]);
yy = (1 - alpha) * ((1 - beta) * corners[1] + beta * corners[7]) +
alpha * ((1 - beta) * corners[3] + beta * corners[5]);
}
return (xx, yy);
}
(xx, yy) - точка на полу комнаты. В цикле вычисляются координаты точки относительно углов точки: (alpha, beta). Имея координаты углов плана можно вычислить где точка (alpha, beta) на плане окажется в комнате. Эта процедура повторяется сто раз и получается ответ.
Отображение которое переводит точку в комнате в точку на плане и снова в точку в комнате - сжимающее. Это значит что если вы возьмёте две любые точки в комнате, отобразите их, расстояние между образами будет меньше чем расстояние между исходными точками. Это очевидно - план всё-таки куда меньше комнаты. Из математического анализа известно что сжимающее отображение имеет предельную точку и она единственная. Вы можете выбрать любую начальную точку, преобразовывать её раз за разом и, после бесконечного количества преобразований точка окажется в пределе. У меня число повторений не бесконечное, но достаточное. Я прикинул что план всяко не больше половины комнаты. Это значит после каждого отображения точка оказывается в два раза ближе к пределу. В начале цикла точка в комнате, то есть, скажем на расстоянии x от предела. В конце цикла точка окажется на расстоянии x / 2^100 от предела. Это меньше точности double, которая равна x / 2^53. Понятно, что если лист по размеру близок к комнате или даже больше комнаты, метод не работает.
Более правильно делать не фиксированное число итераций, а продолжать отображать точку пока не достигнута нужная точность. Оставлю это в качестве упражнения.
Метод называется "решение системы линейных уравнений методом простых итераций".
То что задача имеет такое "дешёвое" решение не означает что вам не нужно знать что такое математический анализ и линейная алгебра. Это решение не во всех ситуациях применимо, в отличие от других ответов.
Если есть точка на отрезке и нужно узнать какой точке она соответствует на другом отрезке с указанными координатами вершин - не нужно делать то, что вы здесь понаворотили. Находим две точки на противоположенных гранях, строим отрезок параллельный другим 2-м граням и на нем находим искомую точку.
s = input()
s_arr = s.split()
width, height = float(s_arr[0]), float(s_arr[1])
s = input()
s_arr = s.split()
corners = [float(s_i) for s_i in s_arr]
# (0,0), (X,0), (X,Y), (0,Y).
# 'долевые' координаты для всех последующих прямоугольников
koefs = []
for i in range(0, 8, 2 ):
koefs.append(corners[i]/width)
koefs.append(corners[i+1]/height)
def get_new_corners(source: list[float], koefs: list[float]) -> list[float]:
dest = []
# (0,0)'
p1_x = source[0] + (source[2] - source[0]) * koefs[0]
p1_y = source[1] + (source[3] - source[1]) * koefs[0]
p2_x = source[6] + (source[4] - source[6]) * koefs[0]
p2_y = source[7] + (source[5] - source[7]) * koefs[0]
p_x = p1_x + (p2_x - p1_x) * koefs[1]
p_y = p1_y + (p2_y - p1_y) * koefs[1]
dest.append(p_x)
dest.append(p_y)
...
...
так 4-е раза для каждой точки ctrl^c ctrl^v и меняем коэффициенты по возрастанию
return dest
while abs(corners[2] - corners[0]) > 1e-15:
new_arr = get_new_corners(corners, koefs)
corners = new_arr
print(corners[0], corners[1])
результат:
10 5
3.0 2.5 1.0 2.5 1.0 1.5 3.0 1.5
2.5000000000000004 2.083333333333333
Не знаю, может это неспортивно, просто дважды вызвать стандартную функцию решения системы линейных уравнений? Но по мне, так получается весьма короткое решение, без всяких циклов и почти из "первых" принципов.
Правда, я так и не нашёл способа, как эффективно применить ограничение "Переговорка и ее план имеют форму прямоугольника". Похоже, проще представить план в более общем виде аффинного преобразования.
- Матрицу аффинного преобразования "переговорная->план" можно получить решая систему уравнений: RoomT AT == CornerT;
- Неподвижную точку: (A - I) fix == 0 и fix2 == 1;
C++ (неявно, в условии точность входных данных - 4 знака):
#include <cassert>
#include <cfloat>
#include <Eigen/Dense>
#include <Eigen/QR>
using namespace Eigen;
Vector2d fixed_point(Vector2d& rm, Matrix<double, 4, 2>& cr) {
// Находит неподвижную точку аффинного преобразования:
// (0,0, rm[0],0, rm[0],rm[1], 0,rm[1]) -> cr
//
// rm - размеры переговорной
// cr - координаты углов плана переговорной внутри переговорной
Matrix<double, 4, 3> rt{{0., 0., 1.},
{rm[0], 0., 1.},
{rm[0], rm[1], 1.},
{0., rm[1], 1.}};
Matrix<double, 4, 3> ct;
ct.block<4, 2>(0, 0) = cr;
ct.col(2).setOnes();
Matrix3d a; // Матрица аффинного преобразования
a = rt.block<3,3>(0,0).householderQr().solve(ct.block<3,3>(0,0))
.transpose();
Matrix3d ai{a}; // (A - I)*f == 0
ai.topLeftCorner(2,2) -= Matrix2d::Identity(); // f[2] == 1
Vector3d f = ai.householderQr().solve(Vector3d::UnitZ());
assert(ct.transpose().isApprox(a*rt.transpose(), 1e-4) && "Нехороший план");
assert(Vector3d::UnitZ().isApprox(ai*f, 20*DBL_EPSILON)
&& "Нехороший план");
assert(f.isApprox(a*f, 20*DBL_EPSILON) && "Нехороший план");
return f.head<2>();
}
Полный текст с тестами, что б не захламлять.
На Python немного атмосфернее получается:
import numpy as np
def fixed_point(rm, cr):
"""Находит неподвижную точку аффинного преобразования:
(0,0, rm[0],0, rm[0],rm[1], 0,rm[1]) -> cr
rm - размеры переговорной
cr - координаты углов плана переговорной внутри переговорной
"""
rt = np.array([[0., 0., 1.],
[rm[0], 0., 1.],
[rm[0], rm[1], 1.],
[0., rm[1], 1.]])
ct = np.ones((4, 3))
ct[:, 0] = cr[0::2];
ct[:, 1] = cr[1::2];
a = np.linalg.solve(rt[:-1], ct[:-1]).T # Матрица аффинного преобразования
ai = a.copy() # (a - I)*fix == 0
ai[:-1,:-1] -= np.eye(2) # fix[2] == 1
f = np.linalg.solve(ai, [0., 0., 1.])
np.testing.assert_allclose(ct.T, [email protected],
rtol=1e-4, atol=1e-4, err_msg="Нехороший план")
np.testing.assert_allclose([0., 0., 1.], ai@f,
rtol=1e-15, atol=1e-15, err_msg="Нехороший план")
np.testing.assert_allclose(f, a@f,
rtol=1e-15, atol=1e-15, err_msg="Нехороший план")
return f[:2]
Jupyter заметка с тестами, их визуализацией и парой формул, что б не захламлять (да и TeX здесь, увы, нет).