Я пытаюсь создать матрицу из 0 значений, где 1 значение заполняет форму эллипса. Мой эллипс был сгенерирован с помощью minVolEllipse.m ( Link 1 ), который возвращает матрицу уравнения эллипса в «центральной форме» и центр эллипса. Затем я использую фрагмент кода из Ellipse_plot.m (из вышеупомянутой ссылки), чтобы параметризовать вектор по большой / меньшей осям, сгенерировать параметрическое уравнение и сгенерировать матрицу преобразованных координат. Вы можете увидеть их код, чтобы увидеть, как это делается. В результате получается матрица, в которой расположены индексы для точек вдоль эллипса. Он не охватывает все значения вдоль контура эллипса, если я не установил число точек сетки N до смехотворно высокого значения.
Когда я использую команды графика или патча MATLAB, я вижу именно тот результат, который мне нужен. Тем не менее, я хочу, чтобы это было представлено в виде матрицы из 0 значений с 1, где патч «заполняет» пробелы. Очевидно, что MATLAB обладает такой функциональностью, но мне еще предстоит найти код для его выполнения. То, что я ищу, похоже на то, как работает bwfill набора инструментов для обработки изображений ( Ссылка 2 ). bwfill не работает для меня, потому что мой эллипс не является смежным, поэтому функция возвращает матрицу, заполненную полностью 1 значениями.
Надеюсь, я обрисовал проблему достаточно хорошо, если нет, пожалуйста, прокомментируйте, и я могу отредактировать пост, чтобы уточнить.
EDIT:
Я разработал стратегию, использующую двухмерный вектор X из Ellipse_plot.m в качестве входных данных для EllipseDirectFit.m ( Ссылка 3 ). Эта функция возвращает коэффициенты для функции эллипса ax ^ 2 + bxy + cy ^ 2 + dx + dy + f = 0. Используя эти коэффициенты, я вычисляю угол между осью X и большой осью эллипса. Этот угол вместе с центральной и большой / вспомогательной осями передается в ellipseMatrix.m ( Link 4 ), который возвращает заполненную матрицу. К сожалению, матрица, кажется, не вращается от того, что я хочу. Вот часть моего кода:
N = 20; %Number of grid points in ellipse
ellipsepoints = clusterpoints(clusterpoints(:,1)==i,2:3)';
[A,C] = minVolEllipse(ellipsepoints,0.001,N);
%%%%%%%%%%%%%%
%
%Adapted from:
% Ellipse_plot.m
% Nima Moshtagh
% nima@seas.upenn.edu
% University of Pennsylvania
% Feb 1, 2007
% Updated: Feb 3, 2007
%%%%%%%%%%%%%%
%
%
% "singular value decomposition" to extract the orientation and the
% axes of the ellipsoid
%------------------------------------
[U D V] = svd(A);
%
% get the major and minor axes
%------------------------------------
a = 1/sqrt(D(1,1))
b = 1/sqrt(D(2,2))
%theta values
theta = [0:1/N:2*pi+1/N];
%
% Parametric equation of the ellipse
%----------------------------------------
state(1,:) = a*cos(theta);
state(2,:) = b*sin(theta);
%
% Coordinate transform
%----------------------------------------
X = V * state;
X(1,:) = X(1,:) + C(1);
X(2,:) = X(2,:) + C(2);
% Output: Elip_Eq = [a b c d e f]' is the vector of algebraic
% parameters of the fitting ellipse:
Elip_Eq = EllipseDirectFit(X')
% http://mathworld.wolfram.com/Ellipse.html gives the equation for finding the angle theta (teta).
% The coefficients from EllipseDirectFit are rescaled to match what is expected in the wolfram link.
Elip_Eq(2) = Elip_Eq(2)/2;
Elip_Eq(4) = Elip_Eq(4)/2;
Elip_Eq(5) = Elip_Eq(5)/2;
if Elip_Eq(2)==0
if Elip_Eq(1) < Elip_Eq(3)
teta = 0;
else
teta = (1/2)*pi;
endif
else
tetap = (1/2)*acot((Elip_Eq(1)-Elip_Eq(3))/(Elip_Eq(2)));
if Elip_Eq(1) < Elip_Eq(3)
teta = tetap;
else
teta = (pi/2)+tetap;
endif
endif
blank_mask = zeros([height width]);
if teta < 0
teta = pi+teta;
endif
%I may need to switch a and b, depending on which is larger (so that the fist is the major axis)
filled_mask1 = ellipseMatrix(C(2),C(1),b,a,teta,blank_mask,1);
РЕДАКТИРОВАТЬ 2:
В ответ на предложение @BenVoigt я написал решение проблемы для цикла здесь:
N = 20; %Number of grid points in ellipse
ellipsepoints = clusterpoints(clusterpoints(:,1)==i,2:3)';
[A,C] = minVolEllipse(ellipsepoints,0.001,N);
filled_mask = zeros([height width]);
for y=0:1:height
for x=0:1:width
point = ([x;y]-C)'*A*([x;y]-C);
if point < 1
filled_mask(y,x) = 1;
endif
endfor
endfor
Хотя это технически решение проблемы, меня интересует неитеративное решение. Я запускаю этот скрипт на многих больших изображениях, и мне нужно, чтобы он был очень быстрым и параллельным.
РЕДАКТИРОВАТЬ 3:
Спасибо @ математическое.coffee за это решение:
[X,Y] = meshgrid(0:width,0:height);
fill_mask=arrayfun(@(x,y) ([x;y]-C)'*A*([x;y]-C),X,Y) < 1;
Однако я считаю, что есть еще лучший способ сделать это. Вот реализация цикла for, которую я сделал, которая работает быстрее, чем обе вышеуказанные попытки:
ellip_mask = zeros([height width]);
[U D V] = svd(A);
a = 1/sqrt(D(1,1));
b = 1/sqrt(D(2,2));
maxab = ceil(max(a,b));
xstart = round(max(C(1)-maxab,1));
xend = round(min(C(1)+maxab,width));
ystart = round(max(C(2)-maxab,1));
yend = round(min(C(2)+maxab,height));
for y = ystart:1:yend
for x = xstart:1:xend
point = ([x;y]-C)'*A*([x;y]-C);
if point < 1
ellip_mask(y,x) = 1;
endif
endfor
endfor
Есть ли способ достичь этой цели (общий размер изображения по-прежнему [ширина по ширине]) без этого цикла for? Причина, по которой это происходит быстрее, заключается в том, что мне не нужно перебирать все изображение, чтобы определить, находится ли моя точка в эллипсе. Вместо этого я могу просто выполнить итерацию по квадратной области, которая является длиной центра +/- наибольшей основной оси.