Провести окружность методом МНК

nata-G

На плоскости имеются N точек.
Необходимо найти такую окружность (т.е. найти координаты x,y её центра и радиус r чтобы сумма квадратов расстояний от каждой точки до этой окружности была минимальной.
Пытался решать в лоб: находить функционал, брать частные производные по r, x,y - получаются крокодилы. Может, есть идеи получше?

syv7

Esli eto nuzhno zaprogat, to mogu posovetovat vospolzovatsa funkciyami iz biblioteki gsl dlya nahozhdeniya minimumov mnogomernih funkciy.
Imho dovolno udobno i bez lishnego gemora .

Ater

А, может быть, имеет смысл перейти к полярным координатам r и phi? Тогда по r надо аппроксимировать константой. А по phi, наверное вообще не надо. Если твоя зависимость смещена относительно нуля, то проще всего найти центр тяжести систкмы точек и вычесть его.

nata-G

попробую полярные координаты...
: Насчет запрогать - надо следать в дельфи, не знаю, есть ли там такая библиотека.

syv7

Nu vmeste s delphi takaya biblioteka tochno ne idet. Bolee togo, vryadli ee kto-nibud perenosil s c++ v pascal.
Tak chto sovetuyu poiskat drugie biblioteki, navernyaka est. (Nu ili Kalitkina vzbotnut' )

Ater

Найти надо, я так понимаю, радиус и центр окружности? Тогда полярные координаты по-любому надо использовать. Кстати, точки с ошибками даны?

nata-G

точки даны без ошибок
все равно, нужны координаты x,y центра окружности. Как их незаметно впихнуть в полярную систему координат?

stm7543347

Тогда по r надо аппроксимировать константой. А по phi, наверное вообще не надо.
И еще, наверное, надо в вариационной постановке указать, что полюс координат тоже ищется...

Ater

Центр окружности находишь, как центр тяжести системы точек: xc = sum x(i)/N, где N -- число точек. С осью y -- аналогично. Вычитаешь центр тяжести из каждой точки и потом просто аппроксимируешь прямой.

Ater

Наверное, можно и не вариационной... А как центр тяжести... При приличном числе точек штука очень устойчивая.

Evgewkin

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=5557&objectType=file
http://www.math.niu.edu/~rusin/known-math/99/circlefit
В архиве програма на Matlab'е с File Exchange

From: mathematik.tu-darmstadt.de (Peter Spellucci)
Subject: Re: Circle fitting
Date: 29 Apr 1999 11:26:19 GMT
Newsgroups: sci.math.num-analysis
In article <7g7u0b$foo$wanadoo.fr>,
"YvesRunfola" <Yves.wanadoo.fr> writes:
|> I'am looking for the simplest method to find the center x0, Y0 and the
|> raduis R of the circle which fit the N data points Xi, Yi i=1,N. (N>3)
|> I wonder if this least square problem can be solved by solving one or more
|> linear systems ?
yes., of course.
your original problem is
\sum{ \sqrtx_i-x0)^2+(y_i-y0)^2) - r }^2 = min_{x0,y0,r},
a nonlnear least squares problem. if the fit is good, i.e. the points (x_i,y_i)
lie almost on a circle, then
\sqrtx_i-x0)^2+(y_i-y0)^2) almost = r ,
hence
\sqrtx_i-x0)^2+(y_i-y0)^2) + r almost = constant = c (say)
multiplying the sum by "c^2" gives
\sum { (x_i-x0)^2+(y_i-y0)^2-r^2 }^2 = min_{x0,y0,r}
with
\alpha=x0^2+y0^2-r^2
this gives the linear least squares problem
\sum { x_i^2+y_i^2 -2*x0*x_i -2*y0*y_i +\alpha}^2 = min_{x0,y0,\alpha}
solving this gives usually a very good estimate for x0,y0,r.
(often this will be already sufficiently good, although not optimal)
some gauss-newton-steps (i.e. solving
linear least problems with a jacobian depending on the
current guess of x0,y0,r) then will suffice to obtain the optimal solution.
software:
see http://plato.la.as.edu/guide.html
hope this helps
peter

nata-G

почему следует брать центр тяжести в качестве центра окружности?
представим, что почти все точки лежат на вертикальном прямой
тогда центр тяжести также лежит на этой прямой, но на самом деле он должен уходить на бесконечность (и радиус окружности тоже бесконечный).
Поэтому идея с центром тяжести не верна.

Nata70

может, проще доделать в полярных?
если уже вычислены координаты центра окружности (хс, ус привести эту окружность в начало координат (т.е., как предлагалось, данный набор точек (х_i, у_i) заменить на (х_i - хс, у_i - ус и перевести каждую точку из набора (х_i - хс, у_i - ус) в набор (r_i, phi_i) в полярных координатах. после чего вычислить среднее из r_i.

nata-G

самое сложное - найти координаты центра окружности.

Nata70

есть идея, основанная на том, что центр окружности - это точка пересечения центральных перпендикуляров любых двух ее (различных) хорд.

zuzaka

естественно если их найти, задача решается устно.
У меня ничего хорошего не получается.
Зачем вам полярная система - убей не понимаю. В декартовой хотя бы исходные соотношения изящные, только при минимизации переменные не до конца развязываются.
И тут у меня два вопроса. Во-первых, кто придумал, что МНК для окружности разрешается единственным образом.
Во-вторых, зачем вообще МНК? Насколько я понимаю, этот метод используется для полиномиальных уравнений по той только причине, что его легко проводить, а вовсе не из-за каких-то особых качеств получающейся кривой. Вроде, МНК ничуть не лучше, скажем, метода наименьшей суммы модулей. Или минимакса.

nata-G

почему МНК - это первое что пришло в голову
Если есть другие правдоподобные оценки - пишите.

Nata70

Идея насчет серединных перпендикуляров не понравилась?
Каиафе: в полярных просто найти радиус окружности (как среднее из r_i).

gala05

использовать генетические алгоритмы
попробовать градиентный спуск

nata-G

идея срединных перперникуляров не нравится вот чем:
допустим, почти все точки хорошо лежать на какой-то окружности. Но незадача: одна пара срединных перпендикуляров все-таки пересекается на бескончности. Это сильно портит картину - получается, что центр исходной окружности тоже уходит на бескочность, если вычислять его как линейную комбинацию координат остальных точек (например, центр масс точек пересечения срединных перпендикуляров)

Nata70

Просто не бери параллельные хорды. (Серединные перпендикуляры параллельных хорд даже в идеальной окружности не дадут однозначного ответа.)

zuzaka

Если найти координаты центра, то радиус окружности в декартовых ищется ничуть не сложнее. А если не найти, то и полярные не помогут.

nata-G

Просто не бери параллельные хорды.
а брать ли хорды, угол между которыми 0.00001 рад?

Nata70

Проверку на равенство нулю угла между хордами действительно придется заменить на неравенство, сравнивающее этот угол с eps. Eps задай как точность углов из набора phi_i.

sirp

попробуй сделать так:
При фиксированных x,y - центре окружности - найти оптимальный радиус и сумму квадратов отклонений (в общем виде). Сумма квадратов отклонений будет зависеть только от чисел R1,R2,...Rn - расстояния от данных точек до этого фиксированного центра x,y. Может, из вида этой функции можно будет придумать идеи, как искать центр. (Мне лень эту функцию рассчитывать)
Оставить комментарий
Имя или ник:
Комментарий: