手書きモデルで固有値解析

# -*- coding: utf-8 -*-
import appex
import numpy as np
from matplotlib import pyplot as plt
from requests_toolbelt.multipart.encoder import MultipartEncoder

import base64
import json
from requests import Request, Session
from io import BytesIO
from PIL import Image

#PILで開いた画像をbase64形式に変換する
def pil_image_to_base64(pil_image):
    buffered = BytesIO()
    pil_image.save(buffered, format="PNG")
    str_encode_file = base64.b64encode(buffered.getvalue()).decode("utf-8")
    return str_encode_file


def pil_image_to_base642(fp):
    str_encode_file = MultipartEncoder(fields={'encoded_image': (fp.name, fp)})
    return str_encode_file


#PILで開いた画像をCloud Vision APIに投げる
def recognize_image(pil_image):
        str_encode_file = pil_image_to_base64(pil_image)
        str_url = "https://vision.googleapis.com/v1/images:annotate?key="
        str_api_key = "取得したAPIキー"
        str_headers = {'Content-Type': 'application/json'}
        str_json_data = {
            'requests': [
                {
                    'image': {
                        'content': str_encode_file
                    },
                    'features': [
                        {
                            #'type': "TEXT_DETECTION",
                            'type': "DOCUMENT_TEXT_DETECTION",
                            'maxResults': 10
                        }
                    ]
                }
            ]
        }

        obj_session = Session()
        obj_request = Request("POST",
                              str_url + str_api_key,
                              data=json.dumps(str_json_data),
                              headers=str_headers
                              )
        obj_prepped = obj_session.prepare_request(obj_request)
        obj_response = obj_session.send(obj_prepped,
                                        verify=True,
                                        timeout=60
                                        )

        if obj_response.status_code == 200:
            with open('data.json', 'w') as outfile:
                json.dump(obj_response.json(), outfile)
                text = get_fullTextAnnotation(obj_response.text)
            return text

        else:
            return "error"
#返ってきたjsonデータの"fullTextAnnotation"部分のテキストを抽出する
def get_fullTextAnnotation(json_data):
    text_dict = json.loads(json_data)
    try:
        text = text_dict["responses"][0]["fullTextAnnotation"]["text"]
        return text
    except:
        print(None)
        return None

# 手書き画像を共有Extensionで取得する
pil_image = appex.get_image()

# Google Cloud Vision API に画像を投げて取得したデータから重量とバネ剛性を取得する
res = recognize_image(pil_image)
sp = res.splitlines()

m = list()
k = list()

print('◼️Google Cloud Vision API の認識結果')
for v in sp:
    v = v.strip()
    v = v.replace(' ', '')
    v = v.replace(',', '.')
    d = v.split('=')
    print(v, d)
    if d[0] == 'm' or d[0] == 'M':
        m.append(float(d[1]))
    if d[0] == 'k' or d[0] == 'K':
        k.append(float(d[1]))

print("◼️バネマス情報")
for v in m:
    print('m =',v)
for v in k:
    print('k =',v)

M1 = len(m)
M2 = M1 + 1
M = np.zeros((M2, M2))
K = np.zeros((M2, M2))

# 質量マトリクスと剛性マトリクスを作成
for i in range(len(m)):
    M[i][i] = m[i] / 9.80665

for i in range(len(k)):
    i1 = i
    i2 = i + 1
    K[i1][i1] = K[i1][i1] + k[i]
    K[i2][i2] = K[i2][i2] + k[i]
    K[i1][i2] = K[i1][i2] - k[i]
    K[i2][i1] = K[i2][i1] - k[i]

# 最後のばねの端部が固定点として行と列を削除して拘束条件とする
M = np.delete(M, M1, 0)
M = np.delete(M, M1, 1)
K = np.delete(K, M1, 0)
K = np.delete(K, M1, 1)

# 質量マトリクスの逆行列を計算
M_inv = np.linalg.inv(M)

# 固有値と固有ベクトルを計算
omega, v = np.linalg.eig(np.dot(M_inv, K))

# 固有値の順番を昇順にソート
freq = list()
omega_sort = np.sort(omega)
for o in omega_sort:
    freq.append(np.sqrt(o)/(2.0*np.pi))

# 固有値のソート時のインデックスを取得
# ⇒固有ベクトルと対応させるため
sort_index = np.argsort(omega)

# 固有値に対応する固有ベクトルをソート
v_sort = []
for i in range(len(sort_index)):
    v_sort.append(v.T[sort_index[i]])
v_sort = np.array(v_sort)

print('◼️固有値解析結果')
# 結果をコンソールに表示
for i in range(len(omega_sort)):
    s = '{0:d}次 {1:10.3f} Hz'.format(i+1,freq[i])
    print(s)
    for j in range(0,v_sort.shape[1]):
        s = '{:10.3f}'.format(v_sort[i][j])
        print(s)

# グラフ化のために自由度軸を作成
dof = np.linspace(len(sort_index), 0, len(sort_index)+1)

# ここからグラフ描画
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 9
plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。
fig = plt.figure(figsize=(3,6))
ax1 = fig.add_subplot(111)
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。
ax1.set_xlabel('Eigen Vector')
ax1.set_ylabel('Degree of freedom')

# データの範囲と刻み目盛を明示する。
ax1.set_xticks(np.arange(-2, 2, 0.5))
ax1.set_yticks(np.arange(0, dof[0]+1, 1))
ax1.set_xlim(-1, 1)
ax1.set_ylim(0, dof[0])

# データプロット 固有ベクトルの形を見る
for i in range(len(sort_index)):
    eigen_vector = np.concatenate([v_sort[i], [0]])
    s = 'Mode' + str(i+1) + ' {:10.3f}Hz'.format(freq[i])
    ax1.plot(eigen_vector, dof, lw=1, marker='o', label=s)

fig.tight_layout()
plt.legend(bbox_to_anchor=(0, 0.1), loc='lower left')

# グラフを表示する。
plt.show()
plt.close()