How to integrate a 2x2 binary matrix in polar coordinates?

I have a matrix of 0s and 1s [512,512] where the 1s form a square in the middle. This matrix is just a test for more complex shapes in the future. What I want to do is compute the double integral of the matrix for theta and rho. What functions should I be using in order to do this? Here is what I have right now:

v = zeros(512,512);

x0 = (512/2); y0 = (512/2);

for k=(y0)/2:(512-(y0)/2) %making the square shape for l=(x0)/2:(512-(x0)/2) v(k,l)=1; end end

fun = @(x,y) v; polarfun = @(theta,r) fun(round(r.*sin(theta)),round(r.*cos(theta))); rmax = @(theta) sqrt(sin(theta).^2 + cos(theta).^2); Q = quad2d(polarfun,0,pi/2,0,rmax)

My goal is to eventually use moment analysis to calculate the symmetry of the shape.