在skimage中绘制渐变椭圆

Question

我想用渐变颜色在 skimage 中绘制一个椭圆蒙版。颜色从椭圆内侧开始变化，在椭圆外侧结束。怎么用skimage或者open-cv画出来？

如下图：

Answer 1

简介

让我们从详细描述示例图像开始。

• 这是一个 4 通道图像（RGB + alpha 透明度），但它只使用灰色阴影。
• 图片非常适合绘图，形状周围只有很小的边距。
• 有一个填充的、抗锯齿的、旋转的外椭圆，被透明背景包围。
• 外椭圆填充了旋转椭圆渐变（同心，与椭圆具有相同的旋转），边缘为黑色，线性渐变为白色在中心。
• 外椭圆覆盖有同心、填充、抗锯齿、旋转的内椭圆（同样旋转，两个轴按相同比例缩放）。填充颜​​色为白色。

此外，让：

• a 和 b 是椭圆的半长轴和半短轴
• theta 是椭圆围绕其中心的旋转角度，以弧度为单位（在 0 a 沿 x 轴的点）
• inner_scale 是 inner 与 outer 椭圆的对应轴之间的比率（即 0.5 表示内部缩小 50%）
• h 和 k 是椭圆中心的 (x,y) 坐标

在演示的代码中，我们将使用以下导入：

import cv2
import numpy as np
import math

并将定义我们绘图的参数（我们将计算h和k）设置为：

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

---

步骤#1

为了生成这样的图像，我们需要采取的第一步是计算我们需要的“画布”（我们将要绘制的图像）的大小。为此，我们可以计算旋转外椭圆的边界框，并在其周围添加一些小边距。

我不知道现有的 OpenCV 函数可以有效地执行此操作，但 StackOverflow 节省了一天——已经有一个相关的 question 与一个 answer 链接到一个有用的 article 讨论这个问题。我们可以利用这些资源提出以下 Python 实现：

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)))

注意：我将浮点大小向上取整，因为我们必须覆盖整个像素，并将左上角和右下角作为整数 (x,y) 对返回。

然后我们可以通过以下方式使用该函数：

# 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

---

步骤 #2

第二步是生成透明层。这是一个单通道 8 位图像，其中黑色 (0) 表示完全透明，白色 (255) 表示完全不透明像素。这个任务相当简单，因为我们可以使用cv2.ellipse。

我们可以以RotatedRect 结构（紧贴椭圆的旋转矩形）的形式定义我们的外椭圆。在 Python 中，这表示为包含以下内容的元组：

• 表示旋转矩形中心的元组（x 和 y 坐标）
• 表示旋转矩形大小（宽度和高度）的元组
• 以度为单位的旋转角度

代码如下：

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)

...以及它产生的图像：

---

步骤#3

作为第三步，我们创建一个单通道（灰度或强度）图像，其中包含我们所需的旋转椭圆渐变。但首先，我们如何使用我们的 (a, b)、theta (θ) 和 (h, k) 参数在数学上根据我们图像的笛卡尔 (x, y) 坐标定义旋转椭圆？

这一次，数学 StackExchange 是拯救世界的工具：question 与我们的问题完全匹配，answer 提供了这个有用的等式：

请注意，对于我们从椭圆中心采取的任何方向，左侧在椭圆周长处的计算结果为 1。它在中心为 0，在周边向 1 线性增加，然后再越过它。

由于没有更好的术语，让我们调用右侧的weight。由于它从中心向外缩放得很好，我们可以用它来计算我们想要的梯度。我们的公式在外部为我们提供白色（在浮点图像的情况下为 1.0），在中心为黑色（0.0）。我们想要倒数，所以我们只需从1.0 中减去weight 并将结果裁剪到[0.0, 1.0] 范围内。

让我们从一个简单的纯 Python 实现开始（如手动迭代代表我们图像的 numpy.array 的各个元素）来计算权重。但是，由于我们是懒惰的程序员，我们将使用 Numpy 将计算出的weights 转换为分级图像，使用矢量减法以及numpy.clip。

代码如下：

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)

...以及它产生的图像：

这一切都很好，但由于我们遍历每个像素并在 Python 解释器中进行计算，所以它也是 aaawwwfffuuullllllyyy ssslllooowww .... 可能需要一秒钟，但我们使用的是 Numpy，所以我们肯定可以做到如果我们利用它会更好。这意味着我们可以尽可能地矢量化。

首先，让我们注意到唯一变化的输入是每个给定像素的坐标。这意味着为了向量化我们的算法，我们需要两个数组（与图像大小相同）作为输入，分别保存每个像素的x 和y 坐标。幸运的是，Numpy 为我们提供了生成此类数组的工具——numpy.mgrid。我们可以写

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

生成我们需要的输入数组。然而，让我们观察一下，我们从不直接使用x 和y——而是我们总是用一个常数来抵消它们。让我们通过生成x-h 和y-k 来避免这种偏移操作...

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

我们可以再次预先计算 4 个常量，除此之外，其余的只是向量化 addition、subtraction、multiplication、division 和 powers，它们都是由 Numpy 提供的向量化的操作（即更快）。

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)

与仅使用 Python 的脚本相比，使用此版本的脚本需要大约 30% 的时间。没什么了不起的，但它产生了相同的结果，而且这个任务似乎是你不必经常做的事情，所以对我来说已经足够了。

如果您[读者]知道更快的方法，请将其发布为答案。

现在我们得到了一个浮点图像，其强度范围在 0.0 和 1.0 之间。为了生成我们的结果，我们想要一个 8 位图像，其值在 0 到 255 之间。

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

---

步骤#4

第四步——画内椭圆。这很简单，我们以前做过。我们只需要适当地缩放轴。

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)

这给了我们以下强度图像：

---

步骤#5

第五步——我们快到了。我们所要做的就是将强度和透明度图层组合成一张 BGRA 图像，然后将其保存为 PNG。

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

注意：对红色、绿色和蓝色使用相同的强度只会给我们带来灰色阴影。

当我们保存结果时，我们得到以下图像：

---

结论

鉴于我已经猜测了您用于生成示例图像的参数，我会说我的脚本的结果非常接近。它的运行速度也相当快——如果你想要更好的东西，你可能无法避免关闭裸机（C、C++ 等）。更聪明的方法，或者 GPU 可能会做得更好。值得尝试...

总而言之，这里有一个小演示，证明此代码也适用于其他旋转：

还有我用来写这个的完整脚本：

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))

---

注意：如果您在这篇文章中发现任何愚蠢的错误、令人困惑的术语或任何其他问题，请随时发表建设性评论，或者直接编辑答案以使其更好。我知道有一些方法可以进一步优化这一点——让读者自己来做这个练习（也许提供和补充答案）。