title: mean value in a sphere
body: <p>I'm trying to calculate the mean value of the pixels inside a circle. In the future this needs to be extended to 3D, but for now a 2D sollution would already help me out. </p>
  
  <p>As can be seen in the image, some pixels are entirely inside the circle, but some are only partly inside the circle. The ones partly in the circle also need to contribute only partly to the mean value. The pixels are square. This will simplify the mathematics I hope. </p>
  
  <p>I can calculate the distance from the pixelcorners to the central point, from this you can find the pixels enterly inside and enterly outside. The rest needs correction. But how to find this correction.</p>
  
  <p>[edit] thanks to Heath Raftery the problem is solved! [/edit]</p>
  
  <p><a href="http://www.wolframalpha.com/input/?i=integral(sqrt(r%5E2-x%5E2))" rel="nofollow noreferrer">the integral of a circle with radius r</a></p>
  
  <p>As an example: I want to know the average pixelvalue of pixels in this circle. I know it is 0.3425, since 34.25% of the circle has a value of 1 and the rest is 0.</p>
  
  <p><a href="https://i.stack.imgur.com/gnZoj.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/gnZoj.png" alt="example circle"></a></p>
  
  <p>Function to check what part of a pixel is in the circle:</p>
  
  <pre><code>function [ a ] = incirc( x,y,r )
  %only handles the top right quadrant of a circle
  if x&lt;0||y&lt;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)&lt;r,a=1;return;end %inside circle
  if min(d)&gt;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&lt;1E-5)=NaN; %to find values that are zero
  inters([~((x(1)&lt;inters(1:2))&amp;(inters(1:2)&lt;x(2))),~((y(1)&lt;inters(3:4))&amp;(inters(3:4)&lt;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&amp;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&amp;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&lt;0||a&gt;1
      error('computational error');
  end
  end
  </code></pre>
  
  <p>Script to test the function</p>
  
  <pre><code>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
  </code></pre>
  
  <blockquote>
    <blockquote>
      <p>result = 0.3425</p>
    </blockquote>
  </blockquote>
  
  <p>You can see the algorithm finds the part of the pixel in the circle.</p>
  
  <p><a href="https://i.stack.imgur.com/6l5fm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/6l5fm.png" alt="weights"></a></p>