Аппроксимация дуги окружности
Я уже, наверное, с неделю пытаюсь набрести на более-менее приличный алгоритм аппроксимации дуги окружности вида (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-связанных областей провести внутри каждой области весовой отбор точек, разделив их на несколько групп в зависимости от интенсивности исходного изображения в конкретных точках; затем для каждой группы найти "свой" центр вращения и свести их воедино, согласно весам.
Необходимость аппроксимировать дуги окружностей возникла у меня в связи с недавними техническими ночами: на новой матрице мы измерили координаты центров вращения звездного поля и поворотного стола (+ просканировали небо на предмет коррекции координатных поправок для системы управления БТА), получилось расхождение в 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-связанных областей провести внутри каждой области весовой отбор точек, разделив их на несколько групп в зависимости от интенсивности исходного изображения в конкретных точках; затем для каждой группы найти "свой" центр вращения и свести их воедино, согласно весам.