ค่าเฉลี่ยในทรงกลม

ฉันกำลังพยายามคำนวณค่าเฉลี่ยของพิกเซลภายในวงกลม ในอนาคตสิ่งนี้จำเป็นต้องขยายเป็น 3D แต่สำหรับตอนนี้ โซลูชัน 2D จะช่วยฉันได้แล้ว

ดังที่เห็นในภาพ บางพิกเซลอยู่ในวงกลมทั้งหมด แต่บางพิกเซลอยู่ในวงกลมเพียงบางส่วนเท่านั้น ส่วนที่อยู่ในวงกลมก็ต้องมีส่วนทำให้ค่าเฉลี่ยเพียงบางส่วนเท่านั้น พิกเซลเป็นสี่เหลี่ยมจัตุรัส นี่จะทำให้คณิตศาสตร์ที่ฉันหวังง่ายขึ้น

ฉันสามารถคำนวณระยะห่างจากมุมพิกเซลถึงจุดศูนย์กลางได้ จากนี้คุณจะพบพิกเซลที่เข้าไปด้านในและด้านนอก ส่วนที่เหลือต้องการการแก้ไข แต่จะหาการแก้ไขนี้ได้อย่างไร

[แก้ไข] ขอบคุณ Heath Raftery ปัญหาได้รับการแก้ไขแล้ว! [/แก้ไข]

อินทิกรัลของวงกลมที่มีรัศมี r< /ก>

ฟังก์ชั่นตรวจสอบว่าส่วนใดของพิกเซลอยู่ในวงกลม:

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

คุณจะเห็นอัลกอริทึมค้นหาส่วนของพิกเซลในวงกลม