eddy_em: (Default)
[personal profile] eddy_em
Я уже, наверное, с неделю пытаюсь набрести на более-менее приличный алгоритм аппроксимации дуги окружности вида (x²+y²=R²), но только сегодня набрел на отличную книгу, написанную автором уймы статей на подобную тематику — Николаем Черновым. Самой приятной неожиданностью было то, что помимо этой замечательной книги, в которой автор исследует множество алгоритмов аппроксимации прямых и окружностей, на его домашней страничке есть и готовый код к книге.


Необходимость аппроксимировать дуги окружностей возникла у меня в связи с недавними техническими ночами: на новой матрице мы измерили координаты центров вращения звездного поля и поворотного стола (+ просканировали небо на предмет коррекции координатных поправок для системы управления БТА), получилось расхождение в 5.8'' между центрами! Захотелось сравнить, таким ли оно было при наблюдениях около года назад. Но тогда наблюдения производились на другой матрице (FLI, 3000х3000), да и для определения центра вращения поворотной платформы мы делали не набор отдельных экспозиций, а одну большую (в результате чего звезды слились в дуги):

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

Вот, кстати, как сработали мои самопальные алгоритмы по определению координат центра вращения звездного поля из набора снимков:

на рисунке отмечены: «o» - опорные объекты на серии снимков для определения центра вращения поворотного стола; «*» - опорные объекты на серии снимков для определения центра вращения звездного поля; «+» - вычисленные мгновенные центры вращения поворотного стола; «•» - вычисленные мгновенные центры вращения звездного поля.

Итак, для начала я решил посмотреть, насколько же хорошо работают описанные Черновым алгоритмы, ну и заодно определить скорость их работы. Для этого я утащил у него те алгоритмы, которые он назвал наилучшими: метод Levenberg-Marquardt (даже не представляю, как по-русски транскрибировать фамилию второго) итерационного вычисления параметров дуги, а также алгебраические методы Пратта и Таубина (Чернов пишет, что последний - наилучший).

Чтобы уменьшить количество опорных точек (да и размер изображения) в опытах, я проредил исходную картинку 3000х3000 в 10 раз, тупо взяв по одному угловому пикселю в каждом квадрате 10х10. Для уменьшения мусора, возникшего из-за моего прореживания, я отфильтровал результат медианой 3x3. Вот все манипуляции с комментариями:

[i h] = fits_read('P2_0001.fit'); % считываем картинку
bw = im2bw(i,5000/65536); % бинаризуем по уровню интенсивности 5000
ima = bw(1:10:3086, 1:10:3103); % прореживаем
imf=medfilt2(ima); % сглаживаем
[L N]=bwlabel(imf); % ищем 8-связанные области
for i = 1:N; idx=find(L==i); printf("%d: %d\n",i, size(idx)); endfor % смотрим, сколько пикселей в каждой
L(find(L==2))=0; % избавляемся от ложной области
L(find(L==3))=2; % переиндексируем
[X Y] = meshgrid(1:size(ima,2), size(ima,1):-1:1); % создаем массивы координат
idx1 = find(L == 1); % номера точек, соответствующих первой области
idx2 = find(L == 2); % -//- второй области
X1=X(idx1);Y1=Y(idx1);X2=X(idx2);Y2=Y(idx2); % координаты точек областей
% с весами я пока не заморачивался, а надо будет


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

Проверим, что нам даст каждый из методов, а заодно определим, сколько времени они выполняются.

tic; c=LM([X1 Y1], [200 200 50]); toc; c
Elapsed time is 0.00286 seconds.
c =
   153.771   155.237    73.864
xap=mean(X1), yap=mean(Y1), rap=sqrt(std(X1)^2+std(Y1)^2),
xap =  196.91
yap =  199.84
rap =  40.128
tic; c=LM([X1 Y1], [xap yap rap]); toc; c
Elapsed time is 0.002 seconds.
c =
   153.771   155.237    73.864
tic; c=LM([X1 Y1], c); toc; c
Elapsed time is 0.000491 seconds.
c =
   153.771   155.237    73.864

Итерационный метод дает одни и те же результаты, независимо от начального приближения. Но, как это и ожидалось, чем ближе начальное приближение к реальным значениям, тем быстрее выполняется алгоритм.


tic; c=PrattSVD([X1 Y1]); toc; c
Elapsed time is 0.000819 seconds.
c =
   153.776   155.245    73.869
tic; c=TaubinSVD([X1 Y1]); toc; c
Elapsed time is 0.000734 seconds.
c =
   153.776   155.245    73.861

Алгебраические методы отработали аж в 3 раза быстрей итерационного, да еще и выдали почти идентичные результаты. В принципе, с точностью до одной десятой все методы дали один и тот же результат: xc=153.8, yc=155.2, r=73.9. По времени выполнения (и учитывая рекомендации Чернова) выберем алгоритм Таубина.

Вот такая получается красота на графике:

phi=[1:360]*pi/180; xC=c(3)*cos(phi)+c(1); yC=c(3)*sin(phi)+c(1);
plot(X1,Y1,'.',xC,yC); axis square



Теперь выполним расчеты и для второй дуги:

tic; c2=TaubinSVD([X2 Y2]); toc; c2
Elapsed time is 0.00088 seconds.
c2 =
   153.212   154.702    78.395
xC2=c2(3)*cos(phi)+c2(1); yC2=c2(3)*sin(phi)+c2(1);
plot(X1,Y1,'.',xC,yC, X2,Y2, '.', xC2,yC2,c(1),c(2),'+',c2(1),c2(2),'*'); axis square



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

Итак, "опыты на кошках" дали вполне приличные результаты. После обеда буду тестировать на людях реальных данных. Понятно, что с пороговой бинаризацией надо работать аккуратней. Лучше, наверное, отфильтровать отдельно каждую дугу (по наиболее оптимальным порогам). Еще один вариант - после выделения 8-связанных областей провести внутри каждой области весовой отбор точек, разделив их на несколько групп в зависимости от интенсивности исходного изображения в конкретных точках; затем для каждой группы найти "свой" центр вращения и свести их воедино, согласно весам.



This account has disabled anonymous posting.
If you don't have an account you can create one now.
HTML doesn't work in the subject.
More info about formatting

If you are unable to use this captcha for any reason, please contact us by email at support@dreamwidth.org

May 2025

S M T W T F S
    123
45678910
11121314151617
1819202122 2324
25262728293031

Most Popular Tags

Style Credit

Expand Cut Tags

No cut tags
Page generated May. 30th, 2025 05:39 pm
Powered by Dreamwidth Studios