Юрий Красильников писал(а):
В качестве координат предлагаю 57,16316 с.ш., 33,09248 в.д., направление - С.
Я тут немножко поэкспериментировал. Немного уточняю.
Координаты - 57,16322 с.ш., 33,09252 в.д. Направление съемки - С.
Фокусное расстояние теперь 1,95.
Привожу полный текст работающей программы на perl. Она выдает координаты в поперечном Меркаторе 505595,6; 6335557,9, что соответствует вышеуказанным географическим.
Цитата:
# -------------------------------------
$pi = 3.141592653589;
# Ширина снимка
$w=3406;
# Фокусное расстояние
$f0=1.95;
# X-координаты на фото
@ax=(1033,1746,2005,2389);
# Количество элементов массивов
$k=@ax;
# Приводим координаты к диапазону [-0.5;0.5]
for($i=0;$i<$k;$i++)
{
@ax[$i] = @ax[$i]/$w - 0.5;
}
# Координаты на местности (Массиыв должны иметь ту же длину, что и @ax)
@px=(5594.02,5646.08,5666.08,5690.99);
@py=(6011.55,6045.25,6047.71,6029.17);
# Создаем два рабочих массива нужной длины
# Массив азимутов на местности
@arx=($k x 0);
# Массив азимутов по фото
@afx=($k x 0);
# Начальная точка прямой поиска
$ptsx0=5634.84;
$ptsy0=5608.71;
# Конечная точка прямой поиска
$ptsx1=5580.37;
$ptsy1=5538.15;
# Количество точек на прямой, в которых будем сопоставлять азимуты
$kt = 100;
# Сумма квадратов отклонений (заведомо большое число)
$dmin = 10000;
# Основной цикл перебора точек на прямой поиска
for($step=0; $step <= $kt; $step++)
{
# Определяем координаты очередной точки на местности
$ptsx=$ptsx0+$step*($ptsx1-$ptsx0)/$kt;
$ptsy=$ptsy0+$step*($ptsy1-$ptsy0)/$kt;
# Вычисляем азимуты объектов на местности из данной точки
for($i=0;$i<$k;$i++)
{
@arx[$i]=atan2(@py[$i]-$ptsy,$ptsx-@px[$i]);
}
# Вычисляем азимуты объектов на снимке
for($i=0;$i<$k;$i++)
{
# Вычисление значения угла
@afx[$i]=atan2(@ax[$i],$f0);
}
# Определяем поворот, соответствующий минимуму суммы квадратов отклонений
# По сути - определяем точку минимума квадратного трехчлена
# Находим коэффициенты трехчлена a,b,c
$a=0;
$b=0;
$c=0;
for($i=0;$i<$k;$i++)
{
$a++;
$b+=-2*(@afx[$i]-@arx[$i]);
$c+=(@afx[$i]-@arx[$i])*(@afx[$i]-@arx[$i]);
}
# Угол поворота (точка минимума)
$x=-$b/(2*$a);
# Сумма квадратов отклонений (значение трехчлена в этой точке)
$d=$c+$x*($b+$x*$a);
# Если надо - печатаем номер точки, сумму квадратов отклонений,
# среднеквадратичное отклонение (в градусах) и координаты точки для экспорта в Excel
# Тогда надо закомментировать все остальные операторы печати
# print $step, ";", $d, ";", 57.3*sqrt($d/$k), ";", $ptsx, ";", $ptsy, "\n";
# Если сумма квадратов отклонений уменьшилась - нашли лучшее положение и запоминаем
if($d < $dmin)
{
$dmin = $d;
$ptsxs=$ptsx;
$ptsys=$ptsy;
# Печатаем координаты точки, сумму квадратов отклонений и среднеквадратичное отклонение в градусах
printf "%f %f %e %e\n", $ptsxs, $ptsys, $dmin, 57.3*sqrt($dmin/$k);
# Печатаем таблицу углов на фото, на местности (в радианах) и значения отклонений (в градусах)
for($i=0;$i<$k;$i++)
{
printf "%f %f %f\n", @afx[$i], @arx[$i]+$x, 57.3*( @afx[$i]-@arx[$i]-$x);
}
}
}
# В конце печатаем координаты лучшей из найденных точек, сумму квадратов отклонений и среднеквадратичное отклонение в градусах
printf "\n%f %f %e %e\n", $ptsxs, $ptsys, $dmin, 57.3*sqrt($dmin/$k);