Cyclic curves

cyclic_curves.py

from __future__ import annotations

from typing import *

import numpy as np
from matplotlib import pyplot as plt


class Point2D:
    def __init__(self, x: float, y: float):
        self.update(x, y)

    def __str__(self) -> str:
        return f"({self[0]}, {self[1]})"

    def __getitem__(self, index: int) -> float:
        return self.__x if index == 0 else self.__y

    def move(self, dx: float, dy: float):
        self.__x += dx
        self.__y += dy

    def update(self, x: float, y: float):
        self.__x = x
        self.__y = y

    def norm2(self) -> float:
        return self.__x**2 + self.__y**2

    def inner(self, other: Point2D) -> float:
        return self[0] * other[0] + self[1] * other[1]

    def __mul__(self, other: Point2D) -> float:
        if isinstance(other, Point2D):
            return self.inner(other)
        return other * self

    def __rmul__(self, value: float) -> Point2D:
        return self.__class__(value * self[0], value * self[1])

    def __add__(self, other: Point2D) -> Point2D:
        return self.__class__(self[0] + other[0], self[1] + other[1])


class Basis:
    def __init__(self, degree: int):
        self.__degree = degree
        self.coeffs = np.ones(self.npts, dtype="float64")
        for i in range(1, self.npts):
            self.coeffs[i] = 2 * i * self.coeffs[i - 1] / (2 * i + 1)

    @property
    def degree(self) -> int:
        return self.__degree

    @property
    def npts(self) -> int:
        return self.__degree + 1

    def eval_vector(self, index: int, nodes: Tuple[float]) -> Tuple[float]:
        nodes = 2 * np.pi * (np.array(nodes) - 0.5)
        nodes = 1 + np.cos(nodes + 2 * np.pi * index / (1 + 2 * self.degree))
        nodes = nodes**self.degree
        nodes *= self.coeffs[self.degree] / (2**self.degree)
        return nodes

    def eval(self, nodes: Tuple[float]) -> Tuple[Tuple[float]]:
        return tuple(
            self.eval_vector(index, nodes) for index in range(1 + 2 * self.degree)
        )


class Curve:
    def __init__(self, basis, points: Tuple[Point2D]):
        self.basis = basis
        self.points = points

    def eval(self, nodes: Tuple[float]) -> Tuple[Point2D]:
        matrix = self.basis.eval(nodes)
        return np.dot(self.points, matrix)


degree = 7
basis = Basis(degree)
plt.figure()
usample = np.linspace(0, 1, 129)
matrix = basis.eval(usample)
for i, vector in enumerate(matrix):
    label = r"$T_{%d,%d}$" % (i, degree)
    plt.plot(usample, vector, label=label)

plt.legend()

thetas = np.linspace(0, 2 * np.pi, 2 * degree + 1, endpoint=False)
ctrlpoints = [Point2D(np.cos(th), np.sin(th)) for th in thetas]


plt.figure()
# curve = Curve(basis, points)
# # points = curve.eval(usample)
points = np.dot(ctrlpoints, matrix)
xvals = [point[0] for point in points]
yvals = [point[1] for point in points]
plt.plot(xvals, yvals, label="curve")
ctrlpoints.append(ctrlpoints[0])
xvals = [point[0] for point in ctrlpoints]
yvals = [point[1] for point in ctrlpoints]
plt.plot(xvals, yvals, label="ctrlpoints", ls="dotted")
plt.legend()
plt.gca().set_aspect("equal")
plt.show()

cyclic_curves_less_nodes.py

from __future__ import annotations

from typing import *

import numpy as np
from matplotlib import pyplot as plt


class Point2D:
    def __init__(self, x: float, y: float):
        self.update(x, y)

    def __str__(self) -> str:
        return f"({self[0]}, {self[1]})"

    def __getitem__(self, index: int) -> float:
        return self.__x if index == 0 else self.__y

    def move(self, dx: float, dy: float):
        self.__x += dx
        self.__y += dy

    def update(self, x: float, y: float):
        self.__x = x
        self.__y = y

    def norm2(self) -> float:
        return self.__x**2 + self.__y**2

    def inner(self, other: Point2D) -> float:
        return self[0] * other[0] + self[1] * other[1]

    def __mul__(self, other: Point2D) -> float:
        if isinstance(other, Point2D):
            return self.inner(other)
        return other * self

    def __rmul__(self, value: float) -> Point2D:
        return self.__class__(value * self[0], value * self[1])

    def __add__(self, other: Point2D) -> Point2D:
        return self.__class__(self[0] + other[0], self[1] + other[1])


class Basis:
    def __init__(self, degree: int):
        self.__degree = degree
        self.coeffs = [1, 1, 8/9, 4/5, 128/175]

    @property
    def degree(self) -> int:
        return self.__degree

    @property
    def npts(self) -> int:
        return self.__degree + 1

    def eval_vector(self, index: int, nodes: Tuple[float]) -> Tuple[float]:
        nodes = np.array(nodes) - 0.5 + (index / (1 + self.degree))
        nodes = 1 + np.cos(2*np.pi * nodes)
        nodes = nodes**self.degree
        nodes *= self.coeffs[self.degree] / (2**self.degree)
        return nodes

    def eval(self, nodes: Tuple[float]) -> Tuple[Tuple[float]]:
        return tuple(
            self.eval_vector(index, nodes) for index in range(1 + self.degree)
        )


class Curve:
    def __init__(self, basis, points: Tuple[Point2D]):
        self.basis = basis
        self.points = points

    def eval(self, nodes: Tuple[float]) -> Tuple[Point2D]:
        matrix = self.basis.eval(nodes)
        return np.dot(self.points, matrix)


degree = 3
basis = Basis(degree)
plt.figure()
usample = np.linspace(0, 1, 129)
matrix = basis.eval(usample)
for i, vector in enumerate(matrix):
    label = r"$T_{%d,%d}$" % (i, degree)
    plt.plot(usample, vector, label=label)

plt.legend()

thetas = np.linspace(0, 2 * np.pi, degree + 1, endpoint=False)
ctrlpoints = [Point2D(np.cos(th), np.sin(th)) for th in thetas]


plt.figure()
# curve = Curve(basis, points)
# # points = curve.eval(usample)
points = np.dot(ctrlpoints, matrix)
xvals = [point[0] for point in points]
yvals = [point[1] for point in points]
plt.plot(xvals, yvals, label="curve")
ctrlpoints.append(ctrlpoints[0])
xvals = [point[0] for point in ctrlpoints]
yvals = [point[1] for point in ctrlpoints]
plt.plot(xvals, yvals, label="ctrlpoints", ls="dotted")
plt.legend()
plt.gca().set_aspect("equal")
plt.xlim(-1.2, 1.2)
plt.ylim(-1.2, 1.2)
plt.show()