I want to draw an ellipse mask in skimage with gradual change color. The color changes starting from inside ellipse and end at outside ellipse. How to draw it with skimage or open-cv?

Like image below:

Introduction

Let's start by describing the sample image in detail.

It is a 4 channel image (RGB + alpha transparency), but it only uses shades of gray.
The image fits the drawing quite well, there's only a minimal margin around the shape.
There is a filled, anti-aliased, rotated outer ellipse, surrounded by a transparent background.
The outer ellipse is filled with a rotated elliptical gradient (concentric, and with same rotation as the ellipse), which is black on the edge, progressing linearly to white in the center.
The outer ellipse is overlaid with a concentric, filled, anti-aliased, rotated inner ellipse (again same rotation, both axis are scaled in the same proportion). The fill colour is white.

Furthermore, let:

a and b be the semi-major and semi-minor axis of the ellipse
theta be the rotation angle of the ellipse around its center in radians (at 0 a points along the x axis)
inner_scale is the ratio between corresponding axes of the inner to the outer ellipse (i.e. 0.5 means inner is scaled down by 50%)
h and k are the (x,y) coordinates of the ellipse center

In the demonstrated code we will use the following imports:

import cv2
import numpy as np
import math

And set the parameters defining our drawing (we will calculate h and k) as:

a, b = (360.0, 200.0) # Semi-major and semi-minor axis
theta = math.radians(40.0) # Ellipse rotation (radians)
inner_scale = 0.6 # Scale of the inner full-white ellipse

Step #1

The first step we need to take in order to generate such image is to calculate what size of a "canvas" (the image we will be drawing into) we need. To do this, we can calculate the bounding box of the rotated outer ellipse and add some small margin around it.

I'm not aware of an existing OpenCV function to efficiently do this, but StackOverflow saves the day -- there already is a related question with an answer linking to a useful article discussing the problem. We can use those resources to come up with the following Python implementation:

def ellipse_bbox(h, k, a, b, theta):
    ux = a * math.cos(theta)
    uy = a * math.sin(theta)
    vx = b * math.cos(theta + math.pi / 2)
    vy = b * math.sin(theta + math.pi / 2)
    box_halfwidth = np.ceil(math.sqrt(ux**2 + vx**2))
    box_halfheight = np.ceil(math.sqrt(uy**2 + vy**2))
    return ((int(h - box_halfwidth), int(k - box_halfheight))
        , (int(h + box_halfwidth), int(k + box_halfheight)))

NB: I round the floating point sizes up, since we have to cover whole pixels, and return top left and bottom right corner as integer (x,y) pairs.

We can then use the function in the following manner:

# Calculate the image size needed to draw this and center the ellipse
_, (h, k) = ellipse_bbox(0, 0, a, b, theta) # Ellipse center
h += 2 # Add small margin
k += 2 # Add small margin
width, height = (h*2+1, k*2+1) # Canvas size

Step #2

Second step is to generate the transparency layer. This is a single-channel 8bit image where black (0) represents fully transparent, and white (255) represents fully opaque pixels. This task is rather simple, since we can use cv2.ellipse.

We can define our outer ellipse in form a RotatedRect structure (a rotated rectangle tightly fitting the ellipse). In Python this is represented as a tuple containing:

tuple representing the center of the rotated rectangle (x and y coordinates)
tuple representing the size of the rotated rectangle (its width and height)
rotation angle in degrees

Here's the code:

ellipse_outer = ((h,k), (a*2, b*2), math.degrees(theta))

transparency = np.zeros((height, width), np.uint8)
cv2.ellipse(transparency, ellipse_outer, 255, -1, cv2.LINE_AA)

... and the image it produces:

Step #3

As third step, we create a single-channel (grayscale, or intensity) image containing our desired rotated elliptical gradient. But first of all, how do we even mathematically define a rotated ellipse in terms of our image's Cartesian (x, y) coordinates, using our (a, b), theta (θ) , and (h, k) parameters?

This time Mathematics StackExchange is the one to save the day: There's a question exactly matching our problem, and an answer providing this useful equation:

Observe that for any direction we take from the center of the ellipse, the left hand side evaluates to 1 at the perimeter of the ellipse. It's 0 at the center, and linearly increases towards the 1 at the perimeter, and then further past it.

Let's call the right-hand side weight for lack of a better term. Since it scales so nicely from the center outwards, we can use it to calculate our desired gradient. Our formula gives us white (1.0 in case of floating point image) on the outside and black (0.0) in the center. We want the inverse, so we simply subtract weight from 1.0 and clip the result to the range [0.0, 1.0].

Let's begin with a simple Python-only implementation (as in manually iterating over the individual elements of the numpy.array representing our image) to calculate the weights. However, as we're lazy programmers, we'll use Numpy for transforming the calculated weights to the grading image, using vectorized subtraction, along with numpy.clip.

Here's the code:

def make_gradient_v1(width, height, h, k, a, b, theta):
    # Precalculate constants
    st, ct =  math.sin(theta), math.cos(theta)
    aa, bb = a**2, b**2

    weights = np.zeros((height, width), np.float64)    
    for y in range(height):
        for x in range(width):
            weights[y,x] = ((((x-h) * ct + (y-k) * st) ** 2) / aa
                + (((x-h) * st - (y-k) * ct) ** 2) / bb)
            
    return np.clip(1.0 - weights, 0, 1)

... and the image it produces:

This is all nice and dandy, but since we iterate over each pixel and do calculations in the Python interpreter, it's also aaawwwfffuuullllllyyy ssslllooowww.... takes maybe a second, but we're using Numpy, so we can surely do better if we take advantage of it. That means vectorizing whatever we can.

First of all, let's notice that the only input that varies is the coordinates of each given pixel. This means that to vectorize our algorithm, we need as an input two arrays (same size as the image), holding the x and y coordinates of each pixel. Fortunately, Numpy provides us with a tool to generate such arrays -- numpy.mgrid. We can write

y,x = np.mgrid[:height,:width]

to generate the input arrays we need. However, let's observe that we never use x and y directly -- rather we always offset them by a constant. Let's avoid that offset operation by generating x-h and y-k...

y,x = np.mgrid[-k:height-k,-h:width-h]

We can again pre-compute the 4 constants, and beyond that, the rest is just vectorized addition, subtraction, multiplication, division and powers, which are all provided by Numpy as vectorized operations (i.e. much faster).

def make_gradient_v2(width, height, h, k, a, b, theta):
    # Precalculate constants
    st, ct =  math.sin(theta), math.cos(theta)
    aa, bb = a**2, b**2
        
    # Generate (x,y) coordinate arrays
    y,x = np.mgrid[-k:height-k,-h:width-h]
    # Calculate the weight for each pixel
    weights = (((x * ct + y * st) ** 2) / aa) + (((x * st - y * ct) ** 2) / bb)

    return np.clip(1.0 - weights, 0, 1)

A script using this version takes about 30% of the time compared to the Python-only one. Nothing stunning, but it produces the same result and this task seems like something you don't have to do frequently, so it's good enough for me.

If you [the reader] know a faster way, please post it as an answer.

Now we've got a floating point image, where intensities range between 0.0 and 1.0. To generate our result we want to have a 8bit image with values between 0 and 255.

intensity = np.uint8(make_gradient_v2(width, height, h, k, a, b, theta) * 255)

Step #4

Fourth step -- draw the inner ellipse. This is simple, we've done it before. We just have to scale the axes appropriately.

ellipse_inner = ((h,k), (a*2*inner_scale, b*2*inner_scale), math.degrees(theta))

cv2.ellipse(intensity, ellipse_inner, 255, -1, cv2.LINE_AA)

That gives us the following intensity image:

Step #5

Fifth step -- we're almost there. All we have to do is combine the intensity and transparency layers into a BGRA image, and then save it as PNG.

result = cv2.merge([intensity, intensity, intensity, transparency])

NB: Using same intensity for red, green and blue gives us only shades of gray.

When we save the result, we get the following image:

Conclusion

Given that I've guesstimated the parameters you used to produce the sample image, I'd say the result of my script is damn close. It runs reasonably fast too -- if you want anything better, you probably can't avoid going close the the bare metal (C, C++, etc.). A smarter approach, or perhaps a GPU could do even better perhaps. Worth experimenting with...

To conclude, here's a little demonstration that this code also works for other rotations:

And the complete script I used to write this:

import cv2
import numpy as np
import math

# ============================================================================

def ellipse_bbox(h, k, a, b, theta):
    ux = a * math.cos(theta)
    uy = a * math.sin(theta)
    vx = b * math.cos(theta + math.pi / 2)
    vy = b * math.sin(theta + math.pi / 2)
    box_halfwidth = np.ceil(math.sqrt(ux**2 + vx**2))
    box_halfheight = np.ceil(math.sqrt(uy**2 + vy**2))
    return ((int(h - box_halfwidth), int(k - box_halfheight))
        , (int(h + box_halfwidth), int(k + box_halfheight)))

# ----------------------------------------------------------------------------
        
# Rotated elliptical gradient - slow, Python-only approach
def make_gradient_v1(width, height, h, k, a, b, theta):
    # Precalculate constants
    st, ct =  math.sin(theta), math.cos(theta)
    aa, bb = a**2, b**2

    weights = np.zeros((height, width), np.float64)    
    for y in range(height):
        for x in range(width):
            weights[y,x] = ((((x-h) * ct + (y-k) * st) ** 2) / aa
                + (((x-h) * st - (y-k) * ct) ** 2) / bb)
            
    return np.clip(1.0 - weights, 0, 1)

# ----------------------------------------------------------------------------
    
# Rotated elliptical gradient - faster, vectorized numpy approach
def make_gradient_v2(width, height, h, k, a, b, theta):
    # Precalculate constants
    st, ct =  math.sin(theta), math.cos(theta)
    aa, bb = a**2, b**2
        
    # Generate (x,y) coordinate arrays
    y,x = np.mgrid[-k:height-k,-h:width-h]
    # Calculate the weight for each pixel
    weights = (((x * ct + y * st) ** 2) / aa) + (((x * st - y * ct) ** 2) / bb)

    return np.clip(1.0 - weights, 0, 1)

# ============================================================================ 

def draw_image(a, b, theta, inner_scale, save_intermediate=False):
    # Calculate the image size needed to draw this and center the ellipse
    _, (h, k) = ellipse_bbox(0,0,a,b,theta) # Ellipse center
    h += 2 # Add small margin
    k += 2 # Add small margin
    width, height = (h*2+1, k*2+1) # Canvas size

    # Parameters defining the two ellipses for OpenCV (a RotatedRect structure)
    ellipse_outer = ((h,k), (a*2, b*2), math.degrees(theta))
    ellipse_inner = ((h,k), (a*2*inner_scale, b*2*inner_scale), math.degrees(theta))

    # Generate the transparency layer -- the outer ellipse filled and anti-aliased
    transparency = np.zeros((height, width), np.uint8)
    cv2.ellipse(transparency, ellipse_outer, 255, -1, cv2.LINE_AA)
    if save_intermediate:
        cv2.imwrite("eligrad-t.png", transparency) # Save intermediate for demo

    # Generate the gradient and scale it to 8bit grayscale range
    intensity = np.uint8(make_gradient_v1(width, height, h, k, a, b, theta) * 255)
    if save_intermediate:
        cv2.imwrite("eligrad-i1.png", intensity) # Save intermediate for demo
    
    # Draw the inter ellipse filled and anti-aliased
    cv2.ellipse(intensity, ellipse_inner, 255, -1, cv2.LINE_AA)
    if save_intermediate:
        cv2.imwrite("eligrad-i2.png", intensity) # Save intermediate for demo

    # Turn it into a BGRA image
    result = cv2.merge([intensity, intensity, intensity, transparency])
    return result

# ============================================================================ 

a, b = (360.0, 200.0) # Semi-major and semi-minor axis
theta = math.radians(40.0) # Ellipse rotation (radians)
inner_scale = 0.6 # Scale of the inner full-white ellipse
    
cv2.imwrite("eligrad.png", draw_image(a, b, theta, inner_scale, True))

# ============================================================================ 

rows = []
for j in range(0, 4, 1):
    cols = []
    for i in range(0, 90, 10):
        tile = np.zeros((170, 170, 4), np.uint8)
        image = draw_image(80.0, 50.0, math.radians(i + j * 90), 0.6)
        tile[:image.shape[0],:image.shape[1]] = image
        cols.append(tile)
    rows.append(np.hstack(cols))

cv2.imwrite("eligrad-m.png", np.vstack(rows))

Note: If you find any stupid mistakes, confusing terminology, or any other issues with this post, feel free to leave a constructive comment, or just outright edit the answer to make it better. I'm aware that there are ways to further optimize this -- let it be up to the reader to do this excercise (and perhaps provide and complimentary answer).