网格/格子/椭圆区域的矩阵

Posted 2021-04-12

我可能会使用错误的术语，但寻求一些帮助。

我想为位于椭圆形状周边内的网格生成一个x，y值数组。

这里有代码：http://people.sc.fsu.edu/~jburkardt/c_src/ellipse_grid/ellipse_grid.html在Python中完成此任务。

但是，为了我的目的，椭圆已经旋转到一定程度。当前的等式没有考虑到这一点，需要一些帮助来解释这种转换，不确定如何更改代码来做到这一点？

我一直在研究np.meshrid功能，所以如果有更好的方法可以做到这一点，请说。

非常感谢。

答案

给定欧几里德平面中的椭圆，其最一般的形式为形式中的二次曲线

f(x,y) = a x^2 + 2b x y + c y^2 + 2d x + 2f y + g,

人们可以通过计算中心(x0,y0)

((cd-bf)/(b^2-ac), (af-bd)/(b^2-ac)) 

（参见Ellipse on MathWorld的方程式19和20）。长轴a_m的长度可以通过同一页面上的等式21来计算。

现在足以找到圆圈内的所有网格点(x,y)，其中心为(x0,y0)，半径为a_m，

sign(f(x,y)) = sign(f(x0,y0)).

另一答案

要在椭圆内生成晶格点，我们必须知道水平线与椭圆相交的位置。

以角度θ旋转的零中心椭圆方程：

 x = a * Cos(t) * Cos(theta) - b * Sin(t) * Sin(theta)   
 y = a * Cos(t) * Sin(theta) + b * Sin(t) * Cos(theta)

为了简化计算，我们可以引入伪角度Fi和幅度M（给定椭圆的常数）

 Fi = atan2(a * Sin(theta), b * Cos(theta))
 M = Sqrt((a * Sin(theta))^2 + (b * Cos(theta))^2)

所以

 y = M * Sin(Fi) * Cos(t) + M * Cos(Fi) * Sin(t)
 y/M = Sin(Fi) * Cos(t) +  Cos(Fi) * Sin(t)
 y/M = Sin(Fi + t) 

并且对于位置y处的给定水平线的解决方案是

 Fi + t = ArcSin( y / M)
 Fi + t = Pi - ArcSin( y / M)
 t1 = ArcSin( y / M) - Fi        //note two values
 t2 = Pi - ArcSin( y / M) - Fi

在第一个等式中替换t的两个值，并获得给定Y的X值，并生成一个格点序列

要获得顶部和底部坐标，请区分y

y' = M * Cos(Fi + t) = 0
th = Pi/2 - Fi
tl = -Pi/2 - Fi

找到相应的y并将它们用作行的起始和结束Y坐标。

import math

def ellipselattice(cx, cy, a, b, theta):
    res = []
    at = a * math.sin(theta)
    bt = b * math.cos(theta)
    Fi = math.atan2(at, bt)
    M = math.hypot(at, bt)
    ta = math.pi/2 - Fi
    tb = -math.pi/2 - Fi
    y0 = at * math.cos(ta) + bt *math.sin(ta)
    y1 = at * math.cos(tb) + bt *math.sin(tb)
    y0, y1 = math.ceil(cy + min(y0, y1)), math.floor(cy + max(y0, y1))
    for y  in range(y0, y1+1):
        t1 = math.asin(y / M) - Fi
        t2 = math.pi - math.asin(y / M) - Fi
        x1 = a * math.cos(t1) * math.cos(theta) - b* math.sin(t1) * math.sin(theta)
        x2 = a * math.cos(t2) * math.cos(theta) - b* math.sin(t2) * math.sin(theta)
        x1, x2 = math.ceil(cx + min(x1, x2)), math.floor(cx + max(x1, x2))
        line = [(x, y) for x in range(x1, x2 + 1)]
        res.append(line)
    return res

print(ellipselattice(0, 0, 4, 3, math.pi / 4))