Circular mask with ratio on boundaries

Question

I know I must show some effort but I really have no idea how to solve this problem.

I know how to create circular masks: https://stackoverflow.com/a/44874588/2681662

In the example in link above a mask is a boolean array. Means it shows if the pixel is under the circle or not which I would like to call discrete.

However I want to know how much of each pixel is under the circle. So basically I would have a float mask array and in the boundary of the circle I would have a value shows how much of the pixel is under the circle (percentage). Which I would like to call continuous.

Please note, the numbers given are not calculated and are eyeball estimate.

If anyone can point me to right direction I will be grateful.

Answer 1

I wish to thank to all of you. Specially @Yves Daoust you pointed me to right direction. I wanted to share my solution for those who face similar problem:

My solution used intersections. If I could find the intersection area between a rectangle and a circle I might be able to consider a pixel as a rectangle and the circle would be, obviously, the circle.

Intersection between two shapes:

For this You can use shapely:

To create a circle shapely uses points and adds some buffer (radius) to it:

from shapely.geometry import Point
circle = Point(centerx, centery).buffer(r)

to create a rectangle shapely provides box:

from shapely.geometry import box
rect = box(minx=min_x, miny=min_y, maxx=max_x, maxy=max_y)

One can calculate numerous properties of each shape (which technically are polygons) such as area and bounding points.

from shapely.geometry import Point
circle = Point(centerx, centery).buffer(r)
print(circle.area)

One can calculate the intersection of two polygons and it would return a polygon:

from shapely.geometry import Point, box
circle = Point(centerx, centery).buffer(r)
rect = box(minx=min_x, miny=min_y, maxx=max_x, maxy=max_y)
intersection = circle.intersection(rect)

A pixel is a square with sides of 1 unit (pixel). So the intersection area of a pixel and any other shape would result in a value [0, 1] and it's what we were looking for.

Code

Please notice I used ellipses instead of circles since it's inclusive.

My Package:

from __future__ import annotations

from typing import Union, Tuple
from shapely.geometry import Point, Polygon, box
from shapely.affinity import scale, rotate
from matplotlib.patches import Ellipse
import numpy as np


class Pixel:
    def __init__(self, x: int, y: int) -> None:
        """
        Creates a 1x1 box object on the given coordinates
        :param x: int
                x coordinate
        :param y: int
                y coordinate
        """
        self.x = x
        self.y = y
        self.body = self.__generate()

    def __str__(self) -> str:
        return f"Pixel(x={self.x}, y={self.y})"

    def __repr__(self) -> str:
        return self.__str__()

    def __generate(self) -> Polygon:
        """returns a 1x1 box on self.x, self.y"""
        return box(minx=self.x, miny=self.y, maxx=self.x + 1, maxy=self.y + 1)


class EllipticalMask:
    def __init__(self, center: Tuple[Union[float, int], Union[float, int]],
                 a: Union[float, int], b: Union[float, int], angle: Union[float, int] = 0) -> None:
        """
        Creates an ellipse object on the given coordinates and is able to calculate a mask with given pixels.
        :param center: tuple
                (x, y) coordinates
        :param a: float or int
                sami-major axis of ellipse
        :param b:  float or int
                sami-minor axis of ellipse
        :param angle:  float or int
                angle of ellipse (counterclockwise)
        """
        self.center = center
        self.a = a
        self.b = b
        self.angle = angle

        self.body = self.__generate()

    def __generate(self) -> Polygon:
        """Returns an ellipse with given parameters"""
        return rotate(
            scale(
                Point(self.center[1], self.center[0]).buffer(1),
                self.a,
                self.b
            ),
            self.angle
        )

    def __extreme_points(self) -> dict:
        """Finds extreme points which the polygon lying in"""
        x, y = self.body.exterior.coords.xy
        return {
            "X": {
                "MIN": np.floor(min(x)), "MAX": np.ceil(max(x))
            },
            "Y": {
                "MIN": np.floor(min(y)), "MAX": np.ceil(max(y))
            }
        }

    def __intersepter_pixels(self) -> list:
        """Creates a list of pixel objects which ellipse is covering"""
        points = self.__extreme_points()
        return [
            Pixel(x, y)
            for x in np.arange(points["X"]["MIN"], points["X"]["MAX"] + 1).astype(int)
            for y in np.arange(points["Y"]["MIN"], points["Y"]["MAX"] + 1).astype(int)
            if x >= 0 and y >= 0
        ]

    def mask(self, shape: tuple) -> np.ndarray:
        """
        Returns a float mask
        :param shape: tuple
                the shape of the mask as (width, height)
        :return: ndarray
        """
        pixels = self.__intersepter_pixels()
        mask = np.zeros(shape).astype(float)
        for pixel in pixels:
            ratio = pixel.body.intersection(self.body).area
            mask[pixel.x][pixel.y] = ratio

        return mask

    def matplotlib_artist(self) -> Ellipse:
        """
        Returns a matplotlib artist
        :return: Ellipse
        """
        e = Ellipse(xy=(self.center[0] - 0.5, self.center[1] - 0.5), width=2 * self.a, height=2 * self.b,
                    angle=90 - self.angle)
        e.set_facecolor('none')
        e.set_edgecolor("red")
        return e


class CircularMask(EllipticalMask):
    def __init__(self, center: Tuple[Union[float, int], Union[float, int]],
                 r: Union[float, int]) -> None:
        """
        Uses ellipse to create a circle
        :param center:  tuple
                (x, y) coordinates
        :param r: float or int
                radius of circle
        """
        super(CircularMask, self).__init__(center, r, r, 0)

Usage:

from myLib import EllipticalMask
from matplotlib import pyplot as plt

m = EllipticalMask((50, 50), 25, 15, 20)

mask = m.mask((100, 100))
e = m.matplotlib_artist()

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
ax.imshow(mask)
ax.add_artist(e)
plt.show()

Results:

Any feedback is appreciated.

Answer 2

If you are not after top accuracy, an easy measure is to compute the signed distance of the center of the pixel to the circle (which is the distance to the center minus the radius). Then if the distance is larger than 1, consider the pixel fully outside; if less than 0, fully inside. And of course in between, the distance tells you a fraction between 0 and 100%.

A more accurate technique is to compute that distance at the four corners of the pixel in order to find the two edges that are crossed (the sign changes between the endpoints). This allows you to construct a polygon (triangle, quadrilateral or pentagon, exceptionally hexagon when four sides are crossed) and compute its area by the shoelace formula. This is quite manageable.

For an exact result, you need to add the area of the circular segment between the oblique side of the polygon and the circle (https://en.wikipedia.org/wiki/Circular_segment). Note that there are difficult corner cases.

Answer 3

Look online for a modified form of Xiaolin Wu's anti-aliasing line algorithm, which performs anti-aliasing for circle rendering.

While you are not doing any anti-aliasing here, you can easily use the alpha value this algorithm calculates as an ersatz measure of arc coverage over your set of pixels.

Note that you need only consider the subset of this set, over which the circle segment falls.

Answer 4

I suspect that there isn't an easy formula to give you the exact answer to this problem. If you needed to approximate, you could do something like the following:

Look at all for corners of a square:

1. If all four corners are inside the circle, then the square is inside the circle
2. If all four corners are outside the circle, then the square is outside the circle.
3. Otherwise, divide the square into four equal subsquares and recurse

Recurse as deeply as need be to get whatever accuracy you need.

I hope someone has a better solution.