I need to generate circular masks that represent receptive fields in the human eye. The main problem is that I'm not sure how to calculate the percentage of area covered by pixels that are partially outside the circle. For example, consider this simple code for generating a circular mask:
def circle(size=19, r=7, x_offset=0, y_offset=0):
x0 = y0 = size // 2
x0 = x_offset
y0 = y_offset
y, x = np.ogrid[:size, :size]
y = y[::-1]
return ((x - x0)**2 (y-y0)**2)<= r**2
which generates this circle:
As you can see the output is binary and if the circle covers at least 50% of the area, its weight is set to 1 and if it's less than 50% coverage, it's set to 0. Here is an overlay, showing it more clearly:
Another common method of generating circles that appear more smooth is to use a 2D gaussian:
def gaussian2d(size=19, sigma = 3, amplitude=1, x_offset=0, y_offset=0):
x = np.arange(0, size, 1, float)
y = x[:,np.newaxis][::-1]
x0 = y0 = size // 2
x0 = x_offset
y0 = y_offset
return amplitude*np.exp(-(((x-x0)**2 /(2*sigma**2)) ((y-y0)**2 /(2*sigma**2))))
resulting in an output that looks like this:
This solves the problem with the binary output (weights set to 1 or 0), but the problem is that the weights start decreasing even for pixels that are 100% covered by the circle, which is not correct.
ideally, the mask should look something like this (this is hand-drawn):
The pixels that are 100% covered by the circle are set to 1 and the ones that are partially covered, are set to the percentage area covered by the circle, for example:
The efficiency (both in terms of runtime and memory usage) isn't all that important, I just need a way to calculate the percentage of the area of the pixel, that is covered by the circle. Accuracy isn't super important either and a 1% error should be fine.
CodePudding user response：
A box and a circle may overlap or not. If the four corners are inside the circle, then that's 100% of your "percentage of the area" target.
So the first step is calculate the four distances from the corners to the center of the circle.
di = sqrt((xi-xc)^2 (yi-yc)^2)
If all di are greater or equal to the radius R of the circle then box is outside the circle. Similary, if all corners verify di <= R then the cell is totally inside the circle.
The rest of cases have one to three corners inside the circle.
Let's see this "one corner outside" case:
The point is calculate the coordinates of points P,Q.
Because you know A-B is horizontal you can get the xQ coordinate by using the y coordinate (same for A and B)
xQ= xc sqrt(R^2-(yA-yc)^2)
where xc,yc are the coordinates of the center of the circle. Notice that yQ = yA = yB
Be aware you can get two values, depending on the sign of the sqrt you use. Take the one that is xA <= xQ <= xB.
Similary, using xA you can find xP.
In this moment you know all coordinates.
Pay attention to straight P,Q. The area below it is easy, just some triangles. Or the total area of the cell deducing the triangle P,A,Q.
The area between the straight and the circumference is the area of a circular sector (which you can browse the Internet to find it) minus a triangle.
You'll need to find the angle of this sector. For that, I advise to use atan2 instead of atan, because the former gives an angle in (0, 2pi) range. And you must work with the substraction atan2(f1)-atan2(f2) normalized into (0,pi) range.
You can draw other cases (e.g. three corners out) and using those methods above find the area you need.
For example you have a box with two corners out, but don't know which sides the circle cross. You can test a side by this:
- Get a side (let's call it AB). If it's horizontal you know its
y= yA=yBcoordinate. - Find the
xcoord in circle for thisy(two solutions, x,-x) - Test
xA < x xB. If both solutions fail then the circle doesn't cross this side.
If the side is vertical then use its y coord to find the x one.
If a side has both corners out, then you can avoid testing it.
CodePudding user response：
I haven't fully tested the code yet, but at first glance, it seems to work. I'll post it here in case anyone in the future needs something similar to this. It's near impossible to document and comment on every possible scenario so sorry about that, just use the code as is:
def getAngle(x,y):
if x == 0:
if y>0:
return 90
else:
return 270
angle = np.arctan(y/x)*180/np.pi
if x>0:
if y<0:
angle =360
else:
angle =180
return angle
def getAngleRad(x,y):
if x == 0:
if y>0:
return 0.5*np.pi
else:
return 1.5*np.pi
angle = np.arctan(y/x)
if x>0:
if y<0:
angle =2*np.pi
else:
angle =np.pi
return angle
def circularMask(size=3, r=1, x_offset=0, y_offset=0):
x0 = y0 = size /2.0
x0 = x_offset
y0 = y_offset
img = np.zeros((size,size))
for j in range(img.shape[0]):
y = y0-j
for i in range(img.shape[1]):
x = i-x0
d1 = np.sqrt(x**2 y**2)
d2 = np.sqrt((x 1)**2 y**2)
d3 = np.sqrt(x**2 (y-1)**2)
d4 = np.sqrt((x 1)**2 (y-1)**2)
corners_inside = 0
for d in [d1,d2,d3,d4]:
if d<=r:
corners_inside =1
if corners_inside==4:
img[j,i]=1
if corners_inside==3 or corners_inside==1:
if x>0 and y>0:
if corners_inside==3:
xp = x 1
yq = y
else:
xp = x
yq = y-1
xq = np.sqrt(r**2-yq**2)
yp = np.sqrt(r**2-xp**2)
traingle_area = 0.5*(xp-xq)*(yq-yp)
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
if corners_inside==3:
img[j,i]=1-traingle_area circle_segment_area
else:
img[j,i]=traingle_area circle_segment_area
if x<0 and y>0:
if corners_inside==3:
xp = x
yq = y
else:
xp = x 1
yq = y-1
xq = -np.sqrt(r**2-yq**2)
yp = np.sqrt(r**2-xp**2)
traingle_area = 0.5*(xq-xp)*(yq-yp)
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
if corners_inside==3:
img[j,i]=1-traingle_area circle_segment_area
else:
img[j,i]=traingle_area circle_segment_area
if x<0 and y<0:
if corners_inside==3:
xp = x
yq = y-1
else:
xp = x 1
yq = y
xq = -np.sqrt(r**2-yq**2)
yp = -np.sqrt(r**2-xp**2)
traingle_area = 0.5*(xq-xp)*(yp-yq)
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
if corners_inside==3:
img[j,i]=1-traingle_area circle_segment_area
else:
img[j,i]=traingle_area circle_segment_area
if x>0 and y<0:
if corners_inside==3:
xp = x 1
yq = y-1
else:
xp = x
yq = y
xq = np.sqrt(r**2-yq**2)
yp = -np.sqrt(r**2-xp**2)
traingle_area = 0.5*(xp-xq)*(yp-yq)
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
if corners_inside==3:
img[j,i]=1-traingle_area circle_segment_area
else:
img[j,i]=traingle_area circle_segment_area
if corners_inside==2:
angle = getAngle(x,y)
if angle<22.5 or angle>337.5:
yp = y
yq = y-1
xp = np.sqrt(r**2-yp**2)
xq = np.sqrt(r**2-yq**2)
trapezoid_area = (0.5*(xp xq))-x
theta_p = np.arctan(yp/xp)
theta_q = np.arctan(yq/xq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
if 22.5<angle<67.5:
if d1<r:
yp = y
yq = y-1
xp = np.sqrt(r**2-yp**2)
xq = np.sqrt(r**2-yq**2)
trapezoid_area = (0.5*(xp xq))-x
else:
xp = x
xq = x 1
yp = np.sqrt(r**2-xp**2)
yq = np.sqrt(r**2-xq**2)
trapezoid_area = (0.5*(yp yq))-y 1
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
if 67.5<angle<112.5:
xp = x
xq = x 1
yp = np.sqrt(r**2-xp**2)
yq = np.sqrt(r**2-xq**2)
trapezoid_area = (0.5*(yp yq))-y 1
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
if 112.5<angle<157.5:
if d3< r:
xp = x 1
xq = x
yp = np.sqrt(r**2-xp**2)
yq = np.sqrt(r**2-xq**2)
trapezoid_area = (0.5*(yp yq))-y 1
else:
yp = y
yq = y-1
xp = -np.sqrt(r**2-yp**2)
xq = -np.sqrt(r**2-yq**2)
trapezoid_area = x 1-(0.5*(xp xq))
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
if 157.5<angle<202.5:
yp = y
yq = y-1
xp = -np.sqrt(r**2-yp**2)
xq = -np.sqrt(r**2-yq**2)
trapezoid_area = x 1-(0.5*(xp xq))
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
if 202.5<angle<247.5:
if d4<r:
yp = y
yq = y-1
xp = -np.sqrt(r**2-yp**2)
xq = -np.sqrt(r**2-yq**2)
trapezoid_area = x 1-(0.5*(xp xq))
else:
xp = x
xq = x 1
yp = -np.sqrt(r**2-xp**2)
yq = -np.sqrt(r**2-xq**2)
trapezoid_area = y-(0.5*(yp yq))
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
if 247.5<angle<292.5:
xp = x
xq = x 1
yp = -np.sqrt(r**2-xp**2)
yq = -np.sqrt(r**2-xq**2)
trapezoid_area = y-(0.5*(yp yq))
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
if 292<angle<337.5:
if d3<r:
yp = y
yq = y-1
xp = np.sqrt(r**2-yp**2)
xq = np.sqrt(r**2-yq**2)
trapezoid_area = (0.5*(xp xq))-x
else:
xq = x
xp = x 1
yp = -np.sqrt(r**2-xp**2)
yq = -np.sqrt(r**2-xq**2)
trapezoid_area = y-(0.5*(yp yq))
theta_p = getAngleRad(xp,yp)
theta_q = getAngleRad(xq,yq)
delta_theta = np.abs(theta_p - theta_q)
circle_segment_area = 0.5*(delta_theta-np.sin(delta_theta))*(r**2)
img[j,i] = trapezoid_area circle_segment_area
return img
the output matches the expected result for the values posted in the question:
Edit: I've found a bug already and will update the code when I manage to fix it.