Я пытаюсь вычислить среднее значение пикселей внутри круга. В будущем это нужно расширить до 3D, но пока решение 2D мне уже помогло бы.
Как видно на изображении, некоторые пиксели полностью находятся внутри круга, а некоторые только частично. Те, что частично в круге, также должны вносить лишь частичный вклад в среднее значение. Пиксели квадратные. Надеюсь, это упростит математику.
Я могу рассчитать расстояние от углов пикселя до центральной точки, из этого вы можете найти пиксели внутри и снаружи. Остальное требует коррекции. Но как найти эту поправку.
[править] благодаря Хиту Рафтери проблема решена! [/редактировать]
интеграл окружности с радиусом r< /а>
Например: я хочу узнать среднее значение пикселей в этом круге. Я знаю, что это 0,3425, так как 34,25% круга имеет значение 1, а остальные 0.
Функция для проверки того, какая часть пикселя находится в круге:
function [ a ] = incirc( x,y,r )
%only handles the top right quadrant of a circle
if x<0||y<0,error('only positive x,y');end
%integral of sqrt(r^2-x^2) dx
F = @(x,r) (1/2)*(x*sqrt(r^2-x^2)+r^2*atan(x/sqrt(r^2-x^2)));
%find corner locations
x=[x-0.5,x+0.5];
y=[y-0.5,y+0.5];
d = sqrt(x.^2+y.^2); %distance to closed and furthest corner
if max(d)<r,a=1;return;end %inside circle
if min(d)>r,a=0;return;end %outside circle
%intersections with edges (r^2 = x^2+y^2)
inters = [sqrt(r^2-y(1)^2),sqrt(r^2-y(2)^2),sqrt(r^2-x(1)^2),sqrt(r^2-x(2)^2)]; %x(1) x(2) y(1) y(2)
%remove imaginary and out of range intersections
inters(imag(inters)~=0)=NaN;
inters(inters<1E-5)=NaN; %to find values that are zero
inters([~((x(1)<inters(1:2))&(inters(1:2)<x(2))),~((y(1)<inters(3:4))&(inters(3:4)<y(2)))])=NaN;
idx = find(~isnan(inters));
if numel(idx)~=2,error('need two intersections of circle with pixel');end
%check area of pixel inside circumference
if all(idx==[1,2]) %2 intersections on y-edge
a=(F(y(2),r)-F(y(1),r)) - x(1); %area
elseif all(idx==[3,4]) %2 intersections on x-edge
a=(F(x(2),r)-F(x(1),r)) - y(1); %area
elseif all(idx==[1,3]) %one intersection on y-edge one on x-edge (left&bottom)
a=(F(inters(1),r)-F(x(1),r))- (y(1)*(inters(1)-x(1)));
elseif all(idx==[2,4]) %one intersection on y-edge one on x-edge (top&right)
a=(inters(2)-x(1))+(F(x(2),r)-F(inters(2),r))-(y(1)*(x(2)-inters(2)));
else
error('geometry')
end
a=real(a);
if a<0||a>1
error('computational error');
end
end
Скрипт для проверки функции
M = ones(100); %data
M(1:50,:)=0;
pos=[50.2,50];
r = 2;
%calculate what the result should be
h=50-pos(2)+0.5;
A=pi*r^2;
wedge = acos(h/r)/pi;
triangle = h*sqrt(r^2-h^2);
res=(A*wedge-triangle)/A
S=0;N=0;
for i = 1:size(M,1)
for j = 1:size(M,2)
x=abs(j-pos(1));
y=abs(i-pos(2));
n=incirc( x,y,r );
M_(i,j)=n;
S = S+M(i,j)*n;
N = N+n;
end
end
result = S/N
результат = 0,3425
Вы можете видеть, что алгоритм находит часть пикселя в круге.
