```python
import numpy as np
import scipy.misc
from scipy.fftpack import dct, idct
import sys
from PIL import Image
import matplotlib
import matplotlib.pyplot as plt
import random
from tqdm._tqdm_notebook import tqdm_notebook
from scipy.fftpack import dct, idct
import seaborn as sns
from skimage.metrics import structural_similarity as ssim
import pandas as pd
import sympy

%matplotlib inline


class ImageLoader:
    def __init__(self, FILE_PATH):
        self.img = np.array(Image.open(FILE_PATH))
        # 行数
        self.row_blocks_count = self.img.shape[0] // 8
        # 列数
        self.col_blocks_count = self.img.shape[1] // 8

    def get_points(self, POINT):
        Row = random.randint(0, len(self.img) - POINT - 1)
        Col = random.randint(0, len(self.img) - 1)
        return self.img[Row : Row + POINT, Col]

    def get_block(self, col, row):
        return self.img[col * 8 : (col + 1) * 8, row * 8 : (row + 1) * 8]


# plt.rcParams['font.family'] ='sans-serif'#使用するフォント
# plt.rcParams["font.sans-serif"] = "Source Han Sans"
plt.rcParams["font.family"] = "Source Han Sans JP"  # 使用するフォント
plt.rcParams["xtick.direction"] = "in"  # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["ytick.direction"] = "in"  # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["xtick.major.width"] = 1.0  # x軸主目盛り線の線幅
plt.rcParams["ytick.major.width"] = 1.0  # y軸主目盛り線の線幅
plt.rcParams["font.size"] = 12  # フォントの大きさ
plt.rcParams["axes.linewidth"] = 1.0  # 軸の線幅edge linewidth。囲みの太さ
matplotlib.font_manager._rebuild()

MONO_DIR_PATH = "../../Mono/"
AIRPLANE = ImageLoader(MONO_DIR_PATH + "airplane512.bmp")
BARBARA = ImageLoader(MONO_DIR_PATH + "barbara512.bmp")
BOAT = ImageLoader(MONO_DIR_PATH + "boat512.bmp")
GOLDHILL = ImageLoader(MONO_DIR_PATH + "goldhill512.bmp")
LENNA = ImageLoader(MONO_DIR_PATH + "lenna512.bmp")
MANDRILL = ImageLoader(MONO_DIR_PATH + "mandrill512.bmp")
MILKDROP = ImageLoader(MONO_DIR_PATH + "milkdrop512.bmp")
SAILBOAT = ImageLoader(MONO_DIR_PATH + "sailboat512.bmp")

IMAGES = [
    AIRPLANE,
    BARBARA,
    BOAT,
    GOLDHILL,
    LENNA,
    MANDRILL,
    MILKDROP,
    SAILBOAT
]
```


```python
class DMLCT:
    def __init__(self, n_bar, N):
        self.n_bar = n_bar
        self.N = N

        self.x_l = (2 * np.arange(N) + 1) / (2 * N)
        self.s_l = np.arange(n_bar) / (n_bar - 1)
        self.xi = (np.arange(n_bar + 1) - 0.5) / (n_bar - 1)

        self.lambda_kh = self.get_lambda_kh(self.n_bar)

        self.w_k_j = self.get_w_k_j(self.n_bar, self.N)
        self.W_L_k_kh = self.get_W_L_k_kh(self.n_bar, self.N)
        self.W_k_kh = self.get_W_k_kh(self.n_bar, self.N)
        self.W_R_k_kh = self.get_W_R_k_kh(self.n_bar, self.N)

    def Lagrange_j(self, j):
        x = sympy.Symbol("x")
        L_x = 1.0
        for l in range(self.n_bar):
            if l != j:
                L_x *= (x - self.s_l[l]) / (self.s_l[j] - self.s_l[l])
        return sympy.integrate(L_x)

    def get_lambda_kh(self, n_bar):
        lambda_kh = np.ones(n_bar)
        lambda_kh[0] = np.sqrt(1 / 2)
        return lambda_kh

    def get_w_k_j(self, n_bar, N):
        L_j = np.zeros((n_bar, N))
        x = sympy.Symbol("x")
        for j in range(n_bar):
            temp = []
            Lj = self.Lagrange_j(j)
            for k in range(N):
                temp.append(Lj.subs(x, self.x_l[k]))
            L_j[j] = np.array(temp)

        w_k_j = np.zeros((n_bar, N))
        for j in range(n_bar):
            w_k_j[j] = scipy.fftpack.dct(L_j[j], norm="ortho")
        return w_k_j

    def get_W_L_k_kh(self, n_bar, N):
        W_L_k_kh = np.zeros((n_bar - 1, N))
        lambda_kh = self.get_lambda_kh(n_bar)

        for kh in range(n_bar - 1):
            W_L_k_kh[kh] = (
                (1 - n_bar)
                * np.sqrt(2 / N)
                * lambda_kh[kh]
                * np.cos(np.pi * kh * (self.xi[0] + 1))
                * self.w_k_j[0]
            )
        return W_L_k_kh

    def get_W_k_kh(self, n_bar, N):
        W_k_kh = np.zeros((n_bar - 1, N))
        for kh in range(n_bar - 1):
            sum_sin = np.zeros(N)
            for j in range(1, n_bar - 2 + 1):
                sum_sin += np.sin(np.pi * kh * self.s_l[j]) * self.w_k_j[j]

            W_k_kh[kh] = (
                (n_bar - 1)
                * np.sqrt(2 / N)
                * self.lambda_kh[kh]
                * (
                    np.cos(np.pi * kh * self.xi[1])
                    * (self.w_k_j[0] - (-1) ** (kh) * self.w_k_j[n_bar - 1])
                    - 2 * np.sin((np.pi * kh) / (2 * (n_bar - 1))) * sum_sin
                )
            )
        return W_k_kh

    def get_W_R_k_kh(self, n_bar, N):
        W_R_k_kh = np.zeros((n_bar - 1, N))
        for kh in range(n_bar - 1):
            W_R_k_kh[kh] = (
                (n_bar - 1)
                * np.sqrt(2 / N)
                * self.lambda_kh[kh]
                * np.cos(np.pi * kh * (self.xi[n_bar] - 1))
                * self.w_k_j[n_bar - 1]
            )
        return W_R_k_kh
```


```python
def get_F_L_k_horizontal(arr, N, row, col):
    # w
    if col == 0:
        w_block = np.zeros(N)
    else:
        w_block = arr[row, (col - 1) * N : col * N]

    return w_block
```


```python
def get_F_R_k_horizontal(arr, N, row, col):
    # e
    if col == arr.shape[1] // N - 1:
        e_block = np.zeros(N)
    else:
        e_block = arr[row, (col + 1) * N : (col + 2) * N]

    return e_block
```


```python
def get_F_L_k_vertical(arr, N, row, col):
    # n
    if row == 0:
        n_block = np.zeros(N)
    else:
        n_block = arr[(row - 1) * N : row * N, col]

    return n_block
```


```python
def get_F_R_k_vertical(arr, N, row, col):
    # s
    if row == arr.shape[0] // N - 1:
        s_block = np.zeros(N)
    else:
        s_block = arr[(row + 1) * N : (row + 2) * N, col]

    return s_block
```


```python
# n_bar = 4
N = 16
```


```python
# dmlct = DMLCT(n_bar, N)
```


```python
IMG = LENNA
```


```python
Fk = np.zeros(IMG.img.shape)
```

# 順変換


```python
def DMLCT_forward(IMG,n_bar,N):
    Fk = np.zeros(IMG.img.shape)
    ## 縦方向
    ### DCT
    for row in range(IMG.img.shape[0] // N):
        for col in range(IMG.img.shape[1]):
            eight_points = IMG.img[N * row : N * (row + 1), col]
            c = scipy.fftpack.dct(eight_points, norm="ortho")
            Fk[N * row : N * (row + 1), col] = c
    ### 残差
    dmlct = DMLCT(n_bar,N)
    for row in range(IMG.img.shape[0] // N):
        for col in range(IMG.img.shape[1]):
            # ビューなら直接いじっちゃう
            F = Fk[N * row : N * (row + 1), col]
            F_L = get_F_L_k_vertical(Fk, N, row, col)
            F_R = get_F_R_k_vertical(Fk, N, row, col)

            U_k_n_bar = np.zeros(N)
            for kh in range(n_bar - 2 + 1):
                U_k_n_bar += (
                    F_L[kh] * dmlct.W_L_k_kh[kh]
                    + F[kh] * dmlct.W_k_kh[kh]
                    + F_R[kh] * dmlct.W_R_k_kh[kh]
                )

            # n_bar = 4 なら 0,1,2は残す 3,4,5,6,7を書き換える
            F[n_bar - 2 + 1 :] -= U_k_n_bar[n_bar - 2 + 1 :]
    # 0を残す
    for k in reversed(range(1, n_bar - 2 + 1)):
        dmlct = DMLCT(k+1, N)
        for row in range(IMG.img.shape[0] // N):
            for col in range(IMG.img.shape[1]):
                # ビューなら直接いじっちゃう
                F = Fk[N * row : N * (row + 1), col]
                F_L = get_F_L_k_vertical(Fk, N, row, col)
                F_R = get_F_R_k_vertical(Fk, N, row, col)

                U_k_n_bar = np.zeros(N)

                for kh in range((k + 1) - 2 + 1):
                    U_k_n_bar += (
                        F_L[kh] * dmlct.W_L_k_kh[kh]
                        + F[kh] * dmlct.W_k_kh[kh]
                        + F_R[kh] * dmlct.W_R_k_kh[kh]
                    )

                F[k] -= U_k_n_bar[k]
    ## 横方向
    ### DCT
    for row in range(Fk.shape[0]):
        for col in range(Fk.shape[1] // N):
            eight_points = Fk[row, N * col : N * (col + 1)]
            c = scipy.fftpack.dct(eight_points, norm="ortho")
            Fk[row, N * col : N * (col + 1)] = c
            
    ### 残差
    dmlct = DMLCT(n_bar,N)
    for row in range(IMG.img.shape[0]):
        for col in range(IMG.img.shape[1] // N):
            F = Fk[row, N * col : N * (col + 1)]
            F_L = get_F_L_k_horizontal(Fk, N, row, col)
            F_R = get_F_R_k_horizontal(Fk, N, row, col)

            U_k_n_bar = np.zeros(N)
            for kh in range(n_bar - 2 + 1):
                U_k_n_bar += (
                    F_L[kh] * dmlct.W_L_k_kh[kh]
                    + F[kh] * dmlct.W_k_kh[kh]
                    + F_R[kh] * dmlct.W_R_k_kh[kh]
                )

            # n_bar = 4 なら 0,1,2は残す 3,4,5,6,7を書き換える
            F[n_bar - 2 + 1 :] -= U_k_n_bar[n_bar - 2 + 1 :]
    # 0を残す
    for k in reversed(range(1, n_bar - 2 + 1)):
        dmlct = DMLCT(k+1, N)
        for row in range(IMG.img.shape[0]):
            for col in range(IMG.img.shape[1] // N):
                F = Fk[row, N * col : N * (col + 1)]
                F_L = get_F_L_k_horizontal(Fk, N, row, col)
                F_R = get_F_R_k_horizontal(Fk, N, row, col)

                U_k_n_bar = np.zeros(N)

                for kh in range((k + 1) - 2 + 1):
                    U_k_n_bar += (
                        F_L[kh] * dmlct.W_L_k_kh[kh]
                        + F[kh] * dmlct.W_k_kh[kh]
                        + F_R[kh] * dmlct.W_R_k_kh[kh]
                    )

                F[k] -= U_k_n_bar[k]
    
    return Fk
```


```python
# DCT係数の平均を求める
```


```python
Fk_values = np.zeros((512,512))
```


```python
for IMG in tqdm_notebook(IMAGES):
    values = np.zeros((25,4))
    Fk = np.zeros(IMG.img.shape)
    for row in range(IMG.img.shape[0] // N):
        for col in range(IMG.img.shape[1] // N):
            block = IMG.img[row * N : (row + 1) * N, col * N : (col + 1) * N]
            c = scipy.fftpack.dct(
                scipy.fftpack.dct(block, axis=0, norm="ortho"), axis=1, norm="ortho"
            )
            Fk[row * N : (row + 1) * N, col * N : (col + 1) * N] = c

    Fk_values += np.abs(Fk)
```


    HBox(children=(IntProgress(value=0, max=8), HTML(value='')))


    



```python
Fk_values /= len(IMAGES)
pd.DataFrame(Fk_values).to_csv("DCT_coef_ave.csv",header=False,index=False)
```


```python
# 各n_barの残差係数の平均を求める
```


```python
Vk_values = np.zeros((512,512))
```


```python
n_bar = 5 
for n_bar in tqdm_notebook(range(2,n_bar+1)):
    Vk_values = np.zeros((512,512))
    for IMG in IMAGES:
        Fk = DMLCT_forward(IMG,n_bar,N)
        Vk_values += np.abs(Fk)
    
    pd.DataFrame(Vk_values / len(IMAGES)).to_csv("DMLCT_" + str(n_bar) + "_coef_ave.csv",header=False,index=False)
```


    HBox(children=(IntProgress(value=0, max=4), HTML(value='')))


    



```python
# DCT係数を読み込む
```


```python
Fk_values = pd.read_csv("DCT_coef_ave.csv",header=None).values
```


```python
# 残差係数を読み込む
```


```python
Vk_values_arr = []
for i in range(2,n_bar+1,1):
    Vk_values_arr.append(pd.read_csv("DMLCT_" + str(i) + "_coef_ave.csv",header=None).values)
```


```python
# NxNブロック1個当たりの係数の平均を求める
```


```python
Fk_block_ave_values = np.zeros((N,N))
for row in range(IMG.img.shape[0] // N):
    for col in range(IMG.img.shape[1] // N):
        if col == 0:
            continue
        if col == IMG.img.shape[1] // N -1:
            continue
        if row == 0:
            continue
        if row == IMG.img.shape[0] // N -1:
            continue
        block = Fk_values[row * N : (row + 1) * N, col * N : (col + 1) * N]
        Fk_block_ave_values += np.abs(block)

Fk_block_ave_values /= (IMG.img.shape[0]//N)**2
```


```python
Vk_block_ave_values_arr = []
for Vk_values in Vk_values_arr:
    Vk_block_ave_values = np.zeros((N,N))
    for row in range(IMG.img.shape[0] // N):
        for col in range(IMG.img.shape[1] // N):
            if col == 0:
                continue
            if col == IMG.img.shape[1] // N -1:
                continue
            if row == 0:
                continue
            if row == IMG.img.shape[0] // N -1:
                continue
            block = Vk_values[row * N : (row + 1) * N, col * N : (col + 1) * N]
            Vk_block_ave_values += np.abs(block)

    Vk_block_ave_values /= (IMG.img.shape[0]//N)**2
    Vk_block_ave_values_arr.append(Vk_block_ave_values)
```


```python
columns = []
for n in reversed(range(2,n_bar+1,1)):
    columns.append("n=" + str(n))
df = pd.DataFrame(columns=columns)
```


```python
k_max = 5
for index in range(1,k_max,1):
    for i in range(index):
        Gk1k2_1 = []
        Gk1k2_2 = []
        for n in reversed(range(2,n_bar+1,1)):
            Vk_block = Vk_block_ave_values_arr[n-2]
            Vk = Vk_block[i,index]
            Fk = Fk_block_ave_values[i,index]
            Gk1k2_1.append(100 * (1 - Vk/Fk))        

            Vk = Vk_block[index,i]
            Fk = Fk_block_ave_values[index,i]
            Gk1k2_2.append(100 * (1 - Vk/Fk))
            
        df.loc["(" + str(i) + "," + str(index) + ")"] = Gk1k2_1
        df.loc["(" + str(index) + "," + str(i) + ")"] = Gk1k2_2
        
    Gk1k2 = []
    for n in reversed(range(2,n_bar+1,1)):
        Vk_block = Vk_block_ave_values_arr[n-2]
        Vk = Vk_block[index,index]
        Fk = Fk_block_ave_values[index,index]
        Gk1k2.append(100 * (1 - Vk/Fk))
    
    df.loc["(" + str(index) + "," + str(index) + ")"] = Gk1k2

Gk1k2 = []
for n in reversed(range(2,n_bar+1,1)):
    Vk_block = Vk_block_ave_values_arr[n-2]
    
    Vk_sum = 0
    Fk_sum = 0
    for row in range(N):
        for col in range(N):
            if row > k_max-1 or col > k_max-1:
                Vk_sum += Vk_block[row,col]
                Fk_sum += Fk_block_ave_values[row,col]
    Gk1k2.append(100 * (1 - Vk_sum/Fk_sum))
df.loc["others"] = Gk1k2
```


```python
df.to_csv("DMLCT_high_freq_comp.csv")
df
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>n=5</th>
      <th>n=4</th>
      <th>n=3</th>
      <th>n=2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>(0,1)</th>
      <td>18.310981</td>
      <td>18.310981</td>
      <td>18.310981</td>
      <td>18.310981</td>
    </tr>
    <tr>
      <th>(1,0)</th>
      <td>21.832817</td>
      <td>21.832817</td>
      <td>21.832817</td>
      <td>21.832817</td>
    </tr>
    <tr>
      <th>(1,1)</th>
      <td>14.575232</td>
      <td>14.575232</td>
      <td>14.575232</td>
      <td>14.575232</td>
    </tr>
    <tr>
      <th>(0,2)</th>
      <td>17.877560</td>
      <td>17.877560</td>
      <td>17.877560</td>
      <td>9.336089</td>
    </tr>
    <tr>
      <th>(2,0)</th>
      <td>16.195630</td>
      <td>16.195630</td>
      <td>16.195630</td>
      <td>8.333530</td>
    </tr>
    <tr>
      <th>(1,2)</th>
      <td>11.987230</td>
      <td>11.987230</td>
      <td>11.987230</td>
      <td>7.455402</td>
    </tr>
    <tr>
      <th>(2,1)</th>
      <td>12.711960</td>
      <td>12.711960</td>
      <td>12.711960</td>
      <td>7.439948</td>
    </tr>
    <tr>
      <th>(2,2)</th>
      <td>12.214650</td>
      <td>12.214650</td>
      <td>12.214650</td>
      <td>5.681808</td>
    </tr>
    <tr>
      <th>(0,3)</th>
      <td>12.405071</td>
      <td>12.405071</td>
      <td>8.628719</td>
      <td>3.997779</td>
    </tr>
    <tr>
      <th>(3,0)</th>
      <td>10.430801</td>
      <td>10.430801</td>
      <td>7.309138</td>
      <td>3.667256</td>
    </tr>
    <tr>
      <th>(1,3)</th>
      <td>7.824086</td>
      <td>7.824086</td>
      <td>5.045261</td>
      <td>2.374177</td>
    </tr>
    <tr>
      <th>(3,1)</th>
      <td>8.051693</td>
      <td>8.051693</td>
      <td>5.461716</td>
      <td>2.826230</td>
    </tr>
    <tr>
      <th>(2,3)</th>
      <td>8.804471</td>
      <td>8.804471</td>
      <td>6.944298</td>
      <td>2.749921</td>
    </tr>
    <tr>
      <th>(3,2)</th>
      <td>7.777447</td>
      <td>7.777447</td>
      <td>5.828321</td>
      <td>2.416075</td>
    </tr>
    <tr>
      <th>(3,3)</th>
      <td>6.360546</td>
      <td>6.360546</td>
      <td>3.840187</td>
      <td>0.963791</td>
    </tr>
    <tr>
      <th>(0,4)</th>
      <td>9.997110</td>
      <td>8.267106</td>
      <td>7.113365</td>
      <td>3.912433</td>
    </tr>
    <tr>
      <th>(4,0)</th>
      <td>7.706722</td>
      <td>6.083833</td>
      <td>5.606183</td>
      <td>3.186106</td>
    </tr>
    <tr>
      <th>(1,4)</th>
      <td>4.039410</td>
      <td>2.487527</td>
      <td>2.614220</td>
      <td>1.213619</td>
    </tr>
    <tr>
      <th>(4,1)</th>
      <td>4.838427</td>
      <td>3.434484</td>
      <td>3.363089</td>
      <td>1.906202</td>
    </tr>
    <tr>
      <th>(2,4)</th>
      <td>5.953920</td>
      <td>4.662677</td>
      <td>4.576186</td>
      <td>2.363430</td>
    </tr>
    <tr>
      <th>(4,2)</th>
      <td>4.229894</td>
      <td>3.112613</td>
      <td>3.123368</td>
      <td>1.878314</td>
    </tr>
    <tr>
      <th>(3,4)</th>
      <td>5.651344</td>
      <td>4.630047</td>
      <td>3.272790</td>
      <td>0.999370</td>
    </tr>
    <tr>
      <th>(4,3)</th>
      <td>4.845592</td>
      <td>3.681487</td>
      <td>2.855581</td>
      <td>0.987264</td>
    </tr>
    <tr>
      <th>(4,4)</th>
      <td>3.622517</td>
      <td>1.985599</td>
      <td>2.198705</td>
      <td>0.982007</td>
    </tr>
    <tr>
      <th>others</th>
      <td>1.054861</td>
      <td>1.299350</td>
      <td>0.865533</td>
      <td>0.328743</td>
    </tr>
  </tbody>
</table>
</div>




```python

```

6be437038f4c5c5dafbb0b9b665db0e6fcde33ae
ipynb
Jupyter Notebook
DMLCT/16x16/DMLCT 交流成分比較_ok.ipynb
Hiroya-W/Python_DCT
5acb7553792335e178d8b99ca1ee42431cc26f92
MIT
DMLCT/16x16/DMLCT 交流成分比較_ok.ipynb
Hiroya-W/Python_DCT
5acb7553792335e178d8b99ca1ee42431cc26f92
MIT
2020-01-06T14:12:55.000Z
2020-02-06T07:00:31.000Z
DMLCT/16x16/DMLCT 交流成分比較_ok.ipynb
Hiroya-W/Python_DCT
5acb7553792335e178d8b99ca1ee42431cc26f92
MIT
Qwen/Qwen-72B
1. YES
2. YES