<reponame>HayatoNatural/Mitsubishi_program
# ここから表面データの処理
import numpy as np
# import matplotlib.pyplot as plt
import pandas as pd
# import os
# import itertools
# from scipy import signal
# import time

class Analysis_surface():

    def __init__(self,k0) -> None:
        self.sampling_num_surface = int(1023) # 表面粗さや2DFFTを計算したりする点数
        if k0 == 1:
            self.resolution = 5.85829e-3 # 表面の測定間隔1, 2DFFTで使用
            self.display_range = 100
        else:
            self.resolution = 0.415229e-3 # 表面の測定間隔2, 2DFFTで使用
            self.display_range = 25
        # self.center = int((self.sampling_num_surface-1)/2) # 表面の中心
        self.k0 = k0
        self.freq_resolution = 1 / (self.resolution*self.sampling_num_surface)
        self.list_i = [1]+list(range(3,22))
        # self.list_k1 = [1,2,3] # ワークの測定場所
        # self.list_k2 = [1,2] # 各場所から何個とるか


    def _read_surface_data(self,i2,k1): # 加工条件i2, 測定場所k1のデータ読み取り(粗い方,1)
        df = pd.read_csv('data/surface_data/'+str(i2)+'-'+str(self.k0)+'-'+str(k1)+'.csv',usecols=range(1023),skipfooter=1,names=list(range(1023)),dtype=float,engine="python") # 生データ
        df2 = df.replace(0,np.nan) # 0を欠損値扱いに
        df3 = df2.interpolate(limit_direction="both") # 欠損値を両側から平均で補完
        z_raw = df3.to_numpy()
        # ここから最小二乗法
        N = self.sampling_num_surface
        x = np.arange(N)*self.resolution
        self.Y,self.X = np.meshgrid(x,x) # 3次元形状でx,y軸を作るときはこれでよい
        X = self.X.reshape(-1,1)
        Y = self.Y.reshape(-1,1)
        z_raw = z_raw.reshape(-1,1)
        ones = np.ones(X.shape)
        Mat = np.hstack([X,Y,ones]) # この行列が最小二乗法の元となる
        # 以下, Ax = b の連立方程式をxについて解く
        A = np.dot(Mat.T,Mat)
        b = np.dot(Mat.T,z_raw)
        x = np.linalg.solve(A,b)
        z_new = z_raw - x[0]*X-x[1]*Y-x[2]
        z_new =z_new.reshape(N,N)
        # z_new = np.mean(df3.to_numpy())
        # df4 = pd.DataFrame(df3.to_numpy() - mean_df3) # 平均面を引いた高さ  Sa = np.mean
        # Sa = np.mean(np.abs(df4.to_numpy()))
        self.surface_data = z_new

    def _caluculate_Sq_Sku(self,surface_data):
        A = ((len(surface_data)-1)*self.resolution)**2
        dA = self.resolution**2
        Sq_2 = np.sum(surface_data**2*dA)/A
        Sq = np.sqrt(Sq_2)
        Sku = np.sum(surface_data**4*dA)/(A*Sq**4)
        return Sq,Sku


    def _do_2DFFT(self):
        n1 = self.display_range
        FFT = np.fft.fft2(self.surface_data) # 変換
        FFT = np.fft.fftshift(FFT) #周波数シフト
        FFT[508:515,508:515] = 1e-3 # 中心に近い低周波成分(±3)を1に
        spec = np.abs(FFT)/(self.sampling_num_surface/2)**2 # パワースペクトル
        spec = spec[511:511+n1,511-n1+1:511+n1] # スペクトルの領域を狭める
        fx = np.arange(n1)*self.freq_resolution
        fy = np.arange(-n1+1,n1)*self.freq_resolution
        FY,FX = np.meshgrid(fy,fx)

        amp = round(np.max(spec),2) # スペクトルの最大値
        idx = np.array(np.unravel_index(np.argmax(spec), spec.shape)) - np.array([0,n1]) #最大値の座標
        print("最大スペクトルの点 : {}".format(idx))
        sp_freq = round(np.sqrt(idx[0]**2+idx[1]**2)*self.freq_resolution,2)  # 最大値の空間周波数(距離に比例)
        fx = round(abs(idx[0])*self.freq_resolution,3)
        fy = round(abs(idx[1])*self.freq_resolution,3)
        angle = round(np.degrees(np.arctan2(fy,fx))) # 最大座標の角度
        return FX,FY,spec,amp,fx,fy,angle


    def _adjust_ax_surface_imshow(self,ax,i2,k1): # 形状データの軸周り
        ax.tick_params(axis = 'x',labelsize =15)
        ax.tick_params(axis = 'y',labelsize =15)
        ax.set_xlabel("X [mm]",size=20,labelpad=5)
        ax.set_ylabel("Y [mm]",size=20,labelpad=5)
        title = "Surface "+str(i2)+"-"+str(self.k0)+"-"+str(k1)
        ax.set_title(title,size=25,y=-0.28)


    def _adjust_ax_2DFFT_imshow(self,ax,i2,k1): # 2DFFTの軸周り
        ax.tick_params(axis = 'x',labelsize =15)
        ax.tick_params(axis = 'y',labelsize =15)
        ax.set_xlabel("FX [1/mm]",size=20,labelpad=5)
        ax.set_ylabel("FY [1/mm]",size=20,labelpad=5)
        # ax.yaxis.set_ticks(np.arange(-8,8))
        title = "2D-FFT "+str(i2)+"-"+str(self.k0)+"-"+str(k1)
        ax.set_title(title,size=25,y=-0.23)

# ここから表面データの処理
import numpy as np
# import matplotlib.pyplot as plt
import pandas as pd
# import os
# import itertools
# from scipy import signal
# import time

class Analysis_surface():

    def __init__(self,k0) -> None:
        self.sampling_num_surface = int(1023) # 表面粗さや2DFFTを計算したりする点数
        if k0 == 1:
            self.resolution = 5.85829e-3 # 表面の測定間隔1, 2DFFTで使用
            self.display_range = 100
        else:
            self.resolution = 0.415229e-3 # 表面の測定間隔2, 2DFFTで使用
            self.display_range = 25
        # self.center = int((self.sampling_num_surface-1)/2) # 表面の中心
        self.k0 = k0
        self.freq_resolution = 1 / (self.resolution*self.sampling_num_surface)
        self.list_i = [1]+list(range(3,22))
        # self.list_k1 = [1,2,3] # ワークの測定場所
        # self.list_k2 = [1,2] # 各場所から何個とるか


    def _read_surface_data(self,i2,k1): # 加工条件i2, 測定場所k1のデータ読み取り(粗い方,1)
        df = pd.read_csv('data/surface_data/'+str(i2)+'-'+str(self.k0)+'-'+str(k1)+'.csv',usecols=range(1023),skipfooter=1,names=list(range(1023)),dtype=float,engine="python") # 生データ
        df2 = df.replace(0,np.nan) # 0を欠損値扱いに
        df3 = df2.interpolate(limit_direction="both") # 欠損値を両側から平均で補完
        z_raw = df3.to_numpy()
        # ここから最小二乗法
        N = self.sampling_num_surface
        x = np.arange(N)*self.resolution
        self.Y,self.X = np.meshgrid(x,x) # 3次元形状でx,y軸を作るときはこれでよい
        X = self.X.reshape(-1,1)
        Y = self.Y.reshape(-1,1)
        z_raw = z_raw.reshape(-1,1)
        ones = np.ones(X.shape)
        Mat = np.hstack([X,Y,ones]) # この行列が最小二乗法の元となる
        # 以下, Ax = b の連立方程式をxについて解く
        A = np.dot(Mat.T,Mat)
        b = np.dot(Mat.T,z_raw)
        x = np.linalg.solve(A,b)
        z_new = z_raw - x[0]*X-x[1]*Y-x[2]
        z_new =z_new.reshape(N,N)
        # z_new = np.mean(df3.to_numpy())
        # df4 = pd.DataFrame(df3.to_numpy() - mean_df3) # 平均面を引いた高さ  Sa = np.mean
        # Sa = np.mean(np.abs(df4.to_numpy()))
        self.surface_data = z_new

    def _caluculate_Sq_Sku(self,surface_data):
        A = ((len(surface_data)-1)*self.resolution)**2
        dA = self.resolution**2
        Sq_2 = np.sum(surface_data**2*dA)/A
        Sq = np.sqrt(Sq_2)
        Sku = np.sum(surface_data**4*dA)/(A*Sq**4)
        return Sq,Sku


    def _do_2DFFT(self):
        n1 = self.display_range
        FFT = np.fft.fft2(self.surface_data) # 変換
        FFT = np.fft.fftshift(FFT) #周波数シフト
        FFT[508:515,508:515] = 1e-3 # 中心に近い低周波成分(±3)を1に
        spec = np.abs(FFT)/(self.sampling_num_surface/2)**2 # パワースペクトル
        spec = spec[511:511+n1,511-n1+1:511+n1] # スペクトルの領域を狭める
        fx = np.arange(n1)*self.freq_resolution
        fy = np.arange(-n1+1,n1)*self.freq_resolution
        FY,FX = np.meshgrid(fy,fx)

        amp = round(np.max(spec),2) # スペクトルの最大値
        idx = np.array(np.unravel_index(np.argmax(spec), spec.shape)) - np.array([0,n1]) #最大値の座標
        print("最大スペクトルの点 : {}".format(idx))
        sp_freq = round(np.sqrt(idx[0]**2+idx[1]**2)*self.freq_resolution,2)  # 最大値の空間周波数(距離に比例)
        fx = round(abs(idx[0])*self.freq_resolution,3)
        fy = round(abs(idx[1])*self.freq_resolution,3)
        angle = round(np.degrees(np.arctan2(fy,fx))) # 最大座標の角度
        return FX,FY,spec,amp,fx,fy,angle


    def _adjust_ax_surface_imshow(self,ax,i2,k1): # 形状データの軸周り
        ax.tick_params(axis = 'x',labelsize =15)
        ax.tick_params(axis = 'y',labelsize =15)
        ax.set_xlabel("X [mm]",size=20,labelpad=5)
        ax.set_ylabel("Y [mm]",size=20,labelpad=5)
        title = "Surface "+str(i2)+"-"+str(self.k0)+"-"+str(k1)
        ax.set_title(title,size=25,y=-0.28)


    def _adjust_ax_2DFFT_imshow(self,ax,i2,k1): # 2DFFTの軸周り
        ax.tick_params(axis = 'x',labelsize =15)
        ax.tick_params(axis = 'y',labelsize =15)
        ax.set_xlabel("FX [1/mm]",size=20,labelpad=5)
        ax.set_ylabel("FY [1/mm]",size=20,labelpad=5)
        # ax.yaxis.set_ticks(np.arange(-8,8))
        title = "2D-FFT "+str(i2)+"-"+str(self.k0)+"-"+str(k1)
        ax.set_title(title,size=25,y=-0.23)
# ここから表面データの処理 # import matplotlib.pyplot as plt # import os # import itertools # from scipy import signal # import time # 表面粗さや2DFFTを計算したりする点数 # 表面の測定間隔1, 2DFFTで使用 # 表面の測定間隔2, 2DFFTで使用 # self.center = int((self.sampling_num_surface-1)/2) # 表面の中心 # self.list_k1 = [1,2,3] # ワークの測定場所 # self.list_k2 = [1,2] # 各場所から何個とるか # 加工条件i2, 測定場所k1のデータ読み取り(粗い方,1) # 生データ # 0を欠損値扱いに # 欠損値を両側から平均で補完 # ここから最小二乗法 # 3次元形状でx,y軸を作るときはこれでよい # この行列が最小二乗法の元となる # 以下, Ax = b の連立方程式をxについて解く # z_new = np.mean(df3.to_numpy()) # df4 = pd.DataFrame(df3.to_numpy() - mean_df3) # 平均面を引いた高さ  Sa = np.mean # Sa = np.mean(np.abs(df4.to_numpy())) # 変換 #周波数シフト # 中心に近い低周波成分(±3)を1に # パワースペクトル # スペクトルの領域を狭める # スペクトルの最大値 #最大値の座標 # 最大値の空間周波数(距離に比例) # 最大座標の角度 # 形状データの軸周り # 2DFFTの軸周り # ax.yaxis.set_ticks(np.arange(-8,8))