本記事では、下図のようなスプライン曲線を作成する雛形コードを載せました。
上図左のスプライン曲線の作成コード部分には、https://teratail.com/questions/279780#reply-398529 を引用しています。本プログラムでは、上図右のようにスプライン曲線上に任意の数の点を生成したり、その点群に対して角度指定で回転させたり、y=x線上に対称に座標変換したデータを生成します。そして、下図のように生成した(x, y)座標データを各々csvファイルに保存できるようにしています（クラスメソッドに追加）。
■本プログラム
import numpy as np import matplotlib.pyplot as plt plt.rcParams['font.size'] = 14 # グラフの基本フォントサイズの設定 import csv class Spline: def __init__(self, xys): # インスタンスの初期化 self.xys = np.array(xys) self.xys = self.xys[self.xys[:, 0].argsort()] # xの昇順にソート N = self.xys.shape[0] assert N >= 3 # 各種パラメータの計算 self.h = np.array([self.xys[i, 0] - self.xys[i - 1, 0] for i in range(1, N)]) self.a = self.xys[:, 1] self.A = np.eye(N) for i in range(2, N): self.A[i - 1, i - 2:i + 1] = self.get_A_line(i) self.r = np.zeros(N) for i in range(2, N): self.r[i - 1] = self.get_r_val(i) def get_A_line(self, i): # 行列Aの各行を計算 return [self.h[i - 2], 2 * (self.h[i - 2] + self.h[i - 1]), self.h[i - 1]] def get_r_val(self, i): # ベクトルrの各要素を計算 return 3 * (self.a[i] - self.a[i - 1]) / self.h[i - 1] - 3 * (self.a[i - 1] - self.a[i - 2]) / self.h[i - 2] def add(self, xy): # 新しい点を追加 self.xys = np.vstack([self.xys, np.array(xy)]) self.xys = self.xys[self.xys[:, 0].argsort()] # 追加後にも再度ソート N = self.xys.shape[0] self.h = np.array([self.xys[i, 0] - self.xys[i - 1, 0] for i in range(1, N)]) self.a = self.xys[:, 1] Anew = np.zeros((N, N)) Anew[0:N - 1, 0:N - 1] = self.A self.A = Anew self.A[N - 2, N - 3:N] = self.get_A_line(N - 1) self.A[N - 1, N - 1] = 1 self.r = np.append(self.r, 0) self.r[N - 2] = self.get_r_val(N - 1) def calc_param(self): # パラメータの計算 Ainv = np.linalg.inv(self.A) self.c = np.dot(Ainv, self.r) N = self.xys.shape[0] self.b = np.array([(self.a[i] - self.a[i - 1]) / self.h[i - 1] - self.h[i - 1] / 3 * (2 * self.c[i - 1] + self.c[i]) for i in range(1, N)]) self.d = np.array([(self.c[i] - self.c[i - 1]) / (3 * self.h[i - 1]) for i in range(1, N)]) def get_y(self, x): # スプライン曲線上のy座標を計算 xs = self.xys[:, 0] if x < xs[0] or x > xs[-1]: return 0 elif x == xs[0]: return self.xys[0, 1] elif x == xs[-1]: return self.xys[-1, 1] i = np.searchsorted(xs, x) - 1 dx = (x - xs[i]) return self.a[i] + self.b[i] * dx + self.c[i] * dx ** 2 + self.d[i] * dx ** 3 def get_points_on_spline(self, num_points): # スプライン曲線上の点を生成 xs = np.linspace(np.min(self.xys[:, 0]), np.max(self.xys[:, 0]), num_points) ys = np.array([self.get_y(x) for x in xs]) return list(zip(xs, ys)) def rotate_points(self, points, angle_degrees): # 点を指定した角度だけ回転させる angle_radians = np.radians(angle_degrees) rotation_matrix = np.array([[np.cos(angle_radians), -np.sin(angle_radians)], [np.sin(angle_radians), np.cos(angle_radians)]]) rotated_points = np.dot(points, rotation_matrix) return rotated_points def reflect_points_y_equals_x(self, points): # y=x線に対して対称な座標変換を行う reflected_points = np.array([[y, x] for x, y in points]) return reflected_points def write_points_to_csv(self, points, filename): # 点の座標をCSVファイルに書き込む with open(filename, 'w', newline='') as csvfile: csv_writer = csv.writer(csvfile) csv_writer.writerow(['X', 'Y']) # ヘッダー行 csv_writer.writerows(points) def main(): ## 乱数の生成 # ランダムな(x, y)座標の生成（x座標を等間隔にし、y座標の上下限を指定） num_points = 10 x_min, x_max = 0, 20 y_min, y_max = 0, 2 if flag_equal_intervals: random_x = np.linspace(x_min, x_max, num_points) # xデータ点を等間隔に生成する else: random_x = np.random.uniform(x_min, x_max, num_points) # x座標を等間隔でなくし、yの上下限を指定 random_y = np.random.uniform(y_min, y_max, num_points) # yの上下限を指定 coordinates = list(zip(random_x, random_y)) print(coordinates) ## ## スプライン曲線の作成 sp = Spline(coordinates) # インスタンスの生成 sp.calc_param() # スプライン曲線の計算 xs = np.linspace(np.min(sp.xys[:, 0]), np.max(sp.xys[:, 0]), 100) ys = np.array([sp.get_y(x) for x in xs]) # スプライン曲線の描画 plt.scatter(sp.xys[:, 0], sp.xys[:, 1], color='blue', s=20, label='raw data') plt.plot(xs, ys, color='red', label='spline') # 作成したスプライン曲線 plt.xlabel('x') plt.ylabel('y') plt.grid() plt.legend() plt.tight_layout() plt.show() # スプライン曲線上の点をいくつか生成 num_points_on_spline = num_points * 5 spline_points = sp.get_points_on_spline(num_points_on_spline) # CSVファイルに出力 csv_filename = "spline_points.csv" sp.write_points_to_csv(spline_points, csv_filename) # スプライン曲線上の点をy=x線に対して対称に変換 reflected_spline_points = sp.reflect_points_y_equals_x(spline_points) # CSVファイルに出力 csv_filename = "reflected_spline_points.csv" sp.write_points_to_csv(reflected_spline_points, csv_filename) # スプライン曲線上の点を回転 rotation_angle_degrees = -90 rotated_spline_points = sp.rotate_points(spline_points, rotation_angle_degrees) # CSVファイルに出力 csv_filename = "rotated_spline_points.csv" sp.write_points_to_csv(rotated_spline_points, csv_filename) plt.scatter(*zip(*spline_points), c='red', s=20, label='spline_points') # スプライン曲線上の点を赤で表示 plt.scatter(*zip(*rotated_spline_points), c='blue', s=20, label='rotated_spline_points') # スプライン曲線上の点を赤で表示 plt.scatter(*zip(*reflected_spline_points), c='magenta', s=20, label='reflected_spline_points') # スプライン曲線上の点を赤で表示 plt.xlabel('x') plt.ylabel('y') plt.grid() plt.legend() plt.tight_layout() plt.show() if __name__ == "__main__": # 乱数シードを固定したい場合は、引数に数値を入れる np.random.seed(99) # x軸に生成する乱数を等間隔にする場合はTrue flag_equal_intervals = True main()
以上
<広告>
リンク
リンク
