Lecture 31. 기저 변환과 영상 압축 — 파이썬 실습

import numpy as np
import plotly.graph_objects as go
from scipy.fft import dct, idct, dctn, idctn

from linalg_viz import COLORS, show_matrix, slider_figure

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(31)
print("numpy", np.__version__)

def 시험사진(n=256, seed=31):
    """skimage 가 없으므로 시험용 이미지를 만든다.

    부드러운 배경 + 또렷한 경계(사각형, 원, 사선) + 잔무늬 + 약한 잡음.
    """
    r = np.random.default_rng(seed)
    y, x = np.mgrid[0:n, 0:n] / (n - 1)
    그림 = 0.34 + 0.24 * (1 - y) * (0.6 + 0.4 * np.cos(1.7 * x))
    그림 += 0.30 * ((x > 0.10) & (x < 0.42) & (y > 0.12) & (y < 0.44))
    그림 -= 0.26 * (np.hypot(x - 0.72, y - 0.32) < 0.16)
    그림 += 0.22 * (np.abs(y - 0.62 * x - 0.30) < 0.020)
    그림 += 0.10 * np.sin(38 * x) * ((y > 0.72) & (y < 0.90)
                                    & (x > 0.15) & (x < 0.85))
    그림 += 0.008 * r.standard_normal((n, n))
    return np.clip(그림, 0, 1)


사진 = 시험사진()
m, n = 사진.shape
오차 = lambda X: float(np.linalg.norm(사진 - X) / np.linalg.norm(사진) * 100)
print(f"{m} x {n} = {사진.size:,} 개의 수")
go.Figure([go.Heatmap(z=사진[::-1], colorscale="gray", zmin=0, zmax=1,
                      showscale=False)],
          layout=dict(title=dict(text="시험 사진"), height=460,
                      xaxis=dict(visible=False, scaleanchor="y"),
                      yaxis=dict(visible=False),
                      margin=dict(l=40, r=40, t=60, b=40)))

def DCT기저(N):
    """열이 DCT-II 기저벡터인 N x N 행렬."""
    자리 = np.arange(N)
    W = np.empty((N, N))
    for k in range(N):
        a = np.sqrt(1 / N) if k == 0 else np.sqrt(2 / N)
        W[:, k] = a * np.cos(np.pi * (자리 + 0.5) * k / N)
    return W

W = DCT기저(8)
print(show_matrix(W, "W  (열이 기저벡터)"))
print("정규직교인가 (W^T W = I) :", np.allclose(W.T @ W, np.eye(8)))
print()
print("0번 기저 (평균) :", np.round(W[:, 0], 4))
print("1번 기저        :", np.round(W[:, 1], 4))
print("각 열의 길이    :", np.round(np.linalg.norm(W, axis=0), 12))
print()
x = rng.normal(size=8)
print("계수는 그냥 내적인가 :", np.allclose(W.T @ x, [W[:, i] @ x for i in range(8)]))
print("scipy 와 같은가      :", np.allclose(W.T @ x, dct(x, norm="ortho")))

x = np.array([7.0, 7.5, 8.0, 8.2, 8.0, 7.0, 5.5, 4.0])
c = W.T @ x
print("x =", x)
print("c =", np.round(c, 4))
print()
print(f"|x| = {np.linalg.norm(x):.9f}")
print(f"|c| = {np.linalg.norm(c):.9f}   <- 파세발")
print("W c 로 되돌아오는가 :", np.allclose(W @ c, x))

def 상위만(v, k):
    """크기가 큰 k 개만 남기고 나머지는 0 으로."""
    남김 = np.zeros_like(v)
    큰것 = np.argsort(-np.abs(v.ravel()))[:k]
    남김.ravel()[큰것] = v.ravel()[큰것]
    return 남김

print(f"{'남길 개수':>10}{'픽셀 기저 오차':>16}{'코사인 기저 오차':>18}"
      f"{'버린 계수 크기':>16}{'같은가':>8}")
for k in (1, 2, 3, 4, 6):
    픽셀 = 상위만(x, k)
    코사인 = 상위만(c, k)
    예측 = float(np.linalg.norm(c - 코사인))
    실제 = float(np.linalg.norm(x - W @ 코사인))
    print(f"{k:>10}{np.linalg.norm(x - 픽셀):>16.6f}{실제:>18.6f}"
          f"{예측:>16.6f}{str(np.isclose(실제, 예측)):>8}")
print()
print("-> 코사인 기저에서는 세 개만 남겨도 오차가 0.59 다. 픽셀 기저는 14.16 이다.")

def 최선인가(B, x, k):
    """큰 것 k 개를 남기는 것이 정말 최선인지 모든 조합과 견준다."""
    직교 = np.allclose(B.T @ B, np.eye(len(x)))
    c = B.T @ x if 직교 else np.linalg.solve(B, x)
    탐욕집합 = frozenset(np.argsort(-np.abs(c))[:k].tolist())
    최선, 최선집합 = np.inf, None
    for 집합 in itertools.combinations(range(len(x)), k):
        남 = np.zeros_like(c)
        남[list(집합)] = c[list(집합)]
        오 = float(np.linalg.norm(x - B @ 남))
        if 오 < 최선:
            최선, 최선집합 = 오, frozenset(집합)
    남 = np.zeros_like(c)
    남[list(탐욕집합)] = c[list(탐욕집합)]
    return float(np.linalg.norm(x - B @ 남)), 최선, 탐욕집합 == 최선집합

import itertools

print("[정규직교 기저] 큰 것부터 고르는 것이 정말 최선인가")
print(f"{'k':>4}{'큰 것 k개':>13}{'모든 조합 중 최선':>18}{'같은가':>8}")
for k in range(1, 6):
    탐, 최, 같 = 최선인가(W, x, k)
    print(f"{k:>4}{탐:>13.6f}{최:>18.6f}{str(같):>8}")
print()
print("-> 늘 최선이다. 서술 파트의 오차 식이 보장해 준다.")

def 기울인기저(W, 세기, rng):
    """정규직교 기저를 무작위로 기울인다. 열 길이는 1 로 유지한다."""
    B = W + 세기 * rng.normal(size=W.shape)
    return B / np.linalg.norm(B, axis=0)

비스듬 = 기울인기저(W, 0.35, np.random.default_rng(1))
print("열 길이가 전부 1 인가 :", np.allclose(np.linalg.norm(비스듬, axis=0), 1))
print("정규직교인가          :", np.allclose(비스듬.T @ 비스듬, np.eye(8)))
print()
print("[비스듬한 기저 하나를 예로]")
print(f"{'k':>4}{'큰 것 k개':>13}{'모든 조합 중 최선':>18}{'같은가':>8}{'몇 배 손해':>12}")
for k in range(1, 6):
    탐, 최, 같 = 최선인가(비스듬, x, k)
    print(f"{k:>4}{탐:>13.6f}{최:>18.6f}{str(같):>8}{탐/최:>12.2f}")

def 실패율(세기, 횟수=200, 씨=7):
    r = np.random.default_rng(씨)
    실패, 전체, 최악 = 0, 0, 1.0
    for _ in range(횟수):
        B = 기울인기저(W, 세기, r)
        if np.linalg.cond(B) > 1e6:
            continue
        for k in range(1, 6):
            탐, 최, 같 = 최선인가(B, x, k)
            전체 += 1
            실패 += (not 같)
            최악 = max(최악, 탐 / max(최, 1e-12))
    return 실패 / 전체, 최악

print(f"{'기울인 정도':>12}{'조건수(중앙값)':>16}{'최선이 아닌 비율':>18}{'최악의 손해':>13}")
r0 = np.random.default_rng(7)
for 세기 in (0.0, 0.05, 0.15, 0.35, 0.60):
    조건 = np.median([np.linalg.cond(기울인기저(W, 세기, r0)) for _ in range(40)])
    비, 최악 = 실패율(세기)
    print(f"{세기:>12.2f}{조건:>16.2f}{비:>18.1%}{최악:>13.2f}")
print()
print("-> 정규직교(0.00)에서는 한 번도 어긋나지 않는다.")
print("   조금만 기울여도 '큰 것부터' 규칙이 무너지기 시작한다.")
print("   비스듬한 기저에서는 어느 것을 버릴지 알려면 다 해 보는 수밖에 없다.")

N = 8
x = rng.normal(size=N)
y = np.concatenate([x, x[::-1]])               # 거울로 이어 붙인다
Y = np.fft.fft(y)
자리 = np.arange(N)

print("거울 확장 :", np.round(y, 3))
print("끝과 시작이 이어지는가 :", np.isclose(y[-1], y[0]), " <- 절벽이 없다")
print()
print(f"{'k':>3}{'Y_k (거울확장의 DFT)':>30}{'유도한 공식':>30}{'같은가':>8}")
for k in range(N):
    공식 = 2 * np.exp(1j*np.pi*k/(2*N)) * np.sum(
        x * np.cos(np.pi * (자리 + 0.5) * k / N))
    print(f"{k:>3}{f'{Y[k].real:+.5f}{Y[k].imag:+.5f}i':>30}"
          f"{f'{공식.real:+.5f}{공식.imag:+.5f}i':>30}"
          f"{str(np.allclose(Y[k], 공식)):>8}")
print()
print("-> DCT 는 새 변환이 아니다. 거울로 이어 붙인 신호의 DFT 다.")

def 감소기울기(신호, 끝=150):
    """계수 크기의 로그-로그 기울기. 진동하는 것은 포락선(국소 최대)으로 잰다."""
    F = np.abs(np.fft.fft(신호))[1:len(신호)//2]
    k = np.arange(1, len(F) + 1)
    봉 = [i for i in range(1, min(끝, len(F) - 1))
          if F[i] > F[i-1] and F[i] > F[i+1] and F[i] > 1e-12]
    if len(봉) < 5:                                  # 진동이 없으면 그냥 잰다
        쓸 = np.flatnonzero(F[:끝] > 1e-12)
        봉 = 쓸.tolist()
    return float(np.polyfit(np.log(k[봉]), np.log(F[봉]), 1)[0])

M = 512
t = np.arange(M)
경우 = (("톱니 (경계에 절벽)", t / M, "절벽 있음"),
        ("계단 하나", (t >= M//3).astype(float), "절벽 있음"),
        ("매끄러운 봉우리", np.exp(-((t - M/2)/40.0)**2), "절벽 없음"))
print(f"{'':>22}{'포락선 기울기':>15}{'':>4}{'뜻':>16}")
for 이름, s, 종류 in 경우:
    기울기 = 감소기울기(s)
    뜻 = "1/k 에 붙는다" if 기울기 > -1.5 else "훨씬 빠르다"
    print(f"{이름:>22}{기울기:>15.3f}{'':>4}{뜻:>16}")
print()
print("-> 절벽이 있으면 -1 근처, 없으면 급전직하. 서술 파트의 부분적분 그대로다.")
print("   계단은 계수가 심하게 진동하므로 그냥 재면 -0.76 이 나온다.")
print("   봉우리만 골라 포락선을 재야 -0.90 으로 -1 에 붙는다.")

램프 = np.linspace(0, 1, 64)
def 필요(계수, 목표=0.999):
    e = np.sort(np.abs(np.asarray(계수)).ravel()**2)[::-1]
    return int(np.searchsorted(np.cumsum(e)/e.sum(), 목표) + 1)
print("곧게 올라가는 기울기 하나를 64점으로 담을 때")
print("  에너지 99.9% 에 필요한 계수")
print("    DFT :", 필요(np.fft.fft(램프)), "개 / 64")
print("    DCT :", 필요(dct(램프, norm='ortho')), "개 / 64")
print()
print("이보다 매끄러운 신호는 없다. 그런데 DFT 는 64개 중 59개를 쓴다.")
print("주기적으로 이어 붙이면서 자기가 만든 절벽을 그리느라 그렇다.")

def 하르행렬(N):
    """서술 파트의 재귀식대로 만든다. 행이 기저벡터."""
    H = np.array([[1.0]])
    while H.shape[0] < N:
        k = H.shape[0]
        H = np.vstack([np.kron(H, [1.0, 1.0]),
                       np.kron(np.eye(k), [1.0, -1.0]) * np.sqrt(k)])
    return H / np.sqrt(N)


def 빠른하르(x):
    """평균과 차이를 반복한다. 덧셈 횟수도 함께 돌려준다."""
    x = np.asarray(x, dtype=float).copy()
    N = 셈 = 0
    N = len(x)
    결과 = np.empty(N)
    자 = np.sqrt(0.5)
    남은 = N
    while 남은 > 1:
        합 = (x[0:남은:2] + x[1:남은:2]) * 자
        차 = (x[0:남은:2] - x[1:남은:2]) * 자
        셈 += 2 * (남은 // 2)
        결과[남은//2:남은] = 차
        x = np.concatenate([합, x[남은:]])
        남은 //= 2
    결과[0] = x[0]
    return 결과, 셈

H8 = 하르행렬(8)
print(show_matrix(H8 * np.sqrt(8), "sqrt(8) H  (행이 기저벡터)"))
print("정규직교인가 :", np.allclose(H8 @ H8.T, np.eye(8)))
print()
print("각 기저벡터가 차지하는 칸 :")
for i, 행 in enumerate(H8):
    자리 = np.flatnonzero(np.abs(행) > 1e-12)
    print(f"  h{i} : {자리.min()}~{자리.max()}  ({len(자리)}칸)")
print()
print("-> 아래로 갈수록 좁아진다. 사인파는 언제나 8칸 전부를 쓴다.")

print("하르 변환의 비용은 O(N) 이다")
print(f"{'N':>7}{'덧셈 횟수':>12}{'2N - 2':>10}{'N log2 N':>12}{'느린 것과 일치':>14}")
for N in (8, 64, 512, 4096):
    v = rng.normal(size=N)
    빠, 셈 = 빠른하르(v)
    느 = 하르행렬(N) @ v
    print(f"{N:>7}{셈:>12,}{2*N-2:>10,}{int(N*np.log2(N)):>12,}"
          f"{str(np.allclose(np.sort(np.abs(빠)), np.sort(np.abs(느)))):>14}")
print()
print("-> FFT 의 N log N 보다도 적다. 조건 ①은 웨이블릿의 압승이다.")

print(f"{'N':>7}{'어느 위치든 최대':>16}{'log2 N + 1':>13}{'DCT 는':>10}")
for N in (64, 256, 1024):
    Hn = 하르행렬(N)
    최대하르 = 최대DCT = 0
    for p in range(1, N, max(1, N // 64)):
        s = np.zeros(N); s[p:] = 1.0
        최대하르 = max(최대하르, int((np.abs(Hn @ s) > 1e-10).sum()))
        최대DCT = max(최대DCT, int((np.abs(dct(s, norm="ortho")) > 1e-10).sum()))
    print(f"{N:>7}{최대하르:>16}{int(np.log2(N))+1:>13}{최대DCT:>10}")
print()
print("-> 하르는 신호가 길어져도 거의 안 늘고, DCT 는 길이만큼 늘어난다.")

# 2차원에서도 마찬가지인가
경계그림 = np.zeros((256, 256)); 경계그림[:, 100:] = 1.0
H = 하르행렬(256)
C경 = dctn(경계그림, norm="ortho")
W경 = H @ 경계그림 @ H.T
print("세로 경계 하나뿐인 256x256 그림 (전체 65,536 개 계수)")
print("  0 이 아닌 DCT  계수 :", int((np.abs(C경) > 1e-10).sum()))
print("  0 이 아닌 하르 계수 :", int((np.abs(W경) > 1e-10).sum()))
print()
# 진짜 사진에서는?
C = dctn(사진, norm="ortho")
Wh = H @ 사진 @ H.T
print("시험 사진에서 같은 개수만 남기면")
print(f"{'남길 비율':>10}{'DCT 오차':>11}{'하르 오차':>12}{'이긴 쪽':>9}")
for 비율 in (0.01, 0.02, 0.05, 0.10):
    k = int(비율 * 사진.size)
    a = 오차(idctn(상위만(C, k), norm="ortho"))
    b = 오차(H.T @ 상위만(Wh, k) @ H)
    print(f"{비율:>10.0%}{a:>11.3f}{b:>12.3f}{('하르' if b<a else 'DCT'):>9}")
print()
print("-> 이 사진은 경계가 많아 하르가 앞선다. JPEG2000 이 웨이블릿을 쓰는 이유다.")

U, s, Vt = np.linalg.svd(사진, full_matrices=False)
SVD랭크 = lambda k: (U[:, :k] * s[:k]) @ Vt[:k]
DCT상위 = lambda k: idctn(상위만(C, k), norm="ortho")

print("[첫 번째 재기] 알맹이 개수 k 를 맞춘다 (사전은 세지 않는다)")
print(f"{'k':>6}{'DCT 상위 k개':>14}{'SVD 랭크 k':>13}{'이긴 쪽':>9}")
for k in (5, 10, 20, 50, 100):
    a, b = 오차(DCT상위(k)), 오차(SVD랭크(k))
    print(f"{k:>6}{a:>14.3f}{b:>13.3f}{('SVD' if b<a else 'DCT'):>9}")
print()
print("-> SVD 의 압승이다. 에카르트-영 정리 그대로다.")

print("[두 번째 재기] 저장할 숫자의 개수를 맞춘다")
print(f"  랭크 1 을 저장하는 값 : U 열 {m}개 + V 열 {n}개 + sigma 1개 = {m+n+1}개")
print(f"  그중 알맹이는 1 개.  사전이 {(m+n)/(m+n+1):.3%}")
print()
print(f"{'예산':>7}{'DCT k':>8}{'SVD 랭크':>9}{'DCT 오차':>10}{'SVD 오차':>10}{'이긴 쪽':>9}")
for 비율 in (0.01, 0.02, 0.05, 0.10, 0.20):
    예산 = int(비율 * 사진.size)
    kd, ks = 예산, max(1, 예산 // (m + n + 1))
    a, b = 오차(DCT상위(kd)), 오차(SVD랭크(ks))
    print(f"{비율:>7.0%}{kd:>8,}{ks:>9}{a:>10.3f}{b:>10.3f}"
          f"{('SVD' if b<a else 'DCT'):>9}")
print()
print("-> 뒤집혔다. 무엇을 예산으로 세느냐가 답을 정한다.")

예산 = int(0.05 * 사진.size)
ks = 예산 // (m + n + 1)
그림 = go.Figure(layout=dict(
    title=dict(text=f"예산 {예산:,}개를 어디에 썼는가"),
    barmode="stack", height=340,
    xaxis=dict(title=dict(text="저장해야 하는 숫자의 개수")),
    margin=dict(l=90, r=30, t=60, b=50)))
그림.add_bar(y=["DCT"], x=[예산], orientation="h", name="계수 (전부 알맹이)",
            marker_color=COLORS["output"])
그림.add_bar(y=["SVD"], x=[ks*m], orientation="h", name="U (사전)",
            marker_color=COLORS["input"])
그림.add_bar(y=["SVD"], x=[ks*n], orientation="h", name="V (사전)",
            marker_color="#2ca02c")
그림.add_bar(y=["SVD"], x=[ks], orientation="h", name="sigma (알맹이)",
            marker_color="#d62728")
그림

print("사진 크기별로 사전이 차지하는 몫")
print(f"{'크기':>16}{'랭크 1 당 저장':>14}{'사전의 몫':>11}")
for (a, b) in ((64, 64), (256, 256), (512, 512), (1920, 1080), (3840, 2160)):
    print(f"{f'{a}x{b}':>16}{a+b+1:>14,}{(a+b)/(a+b+1):>11.4%}")

JPEG표 = np.array([
    [16, 11, 10, 16, 24, 40, 51, 61], [12, 12, 14, 19, 26, 58, 60, 55],
    [14, 13, 16, 24, 40, 57, 69, 56], [14, 17, 22, 29, 51, 87, 80, 62],
    [18, 22, 37, 56, 68, 109, 103, 77], [24, 35, 55, 64, 81, 104, 113, 92],
    [49, 64, 78, 87, 103, 121, 120, 101],
    [72, 92, 95, 98, 112, 100, 103, 99]], float)


def 블록DCT(X, N=8):
    C = np.zeros_like(X)
    for a in range(0, X.shape[0], N):
        for b in range(0, X.shape[1], N):
            C[a:a+N, b:b+N] = dctn(X[a:a+N, b:b+N], norm="ortho")
    return C


def 블록역(C, N=8):
    X = np.zeros_like(C)
    for a in range(0, C.shape[0], N):
        for b in range(0, C.shape[1], N):
            X[a:a+N, b:b+N] = idctn(C[a:a+N, b:b+N], norm="ortho")
    return X


def 양자화(C, q):
    """계수를 간격 q * (JPEG 표) 로 반올림한다."""
    표 = np.tile(JPEG표 * q / 255.0, (C.shape[0]//8, C.shape[1]//8))
    return np.round(C / 표) * 표, 표

블C = 블록DCT(사진)
print("DCT 자체는 손실이 없는가 :", np.allclose(블록역(블C), 사진))
print()
print(f"{'q':>6}{'0 이 된 계수':>13}{'오차':>9}{'경계 낙차 / 안쪽 낙차':>22}")
for q in (0.25, 0.5, 1.0, 2.0, 4.0, 8.0):
    양, 표 = 양자화(블C, q)
    복원 = np.clip(블록역(양), 0, 1)
    죽음 = float((np.abs(양) < 1e-12).mean())
    차 = np.abs(np.diff(복원, axis=1))
    경계, 안쪽 = 차[:, 7::8].mean(), 차[:, [0,1,2,3,4,5]].mean()
    print(f"{q:>6.2f}{죽음:>13.1%}{오차(복원):>9.3f}{경계/안쪽:>22.2f}")
print()
print("-> 세게 걸수록 블록 경계의 낙차가 안쪽보다 커진다. 그것이 격자 자국이다.")
print("   버그가 아니라, 블록을 따로 압축하기로 한 순간 정해진 결과다.")

품질 = [0.15, 0.25, 0.4, 0.6, 1.0, 1.6, 2.5, 4.0, 6.5, 10.0]
조각 = (slice(96, 160), slice(40, 104))
프레임, 이름표 = [], []
for q in 품질:
    양, _ = 양자화(블C, q)
    복원 = np.clip(블록역(양), 0, 1)
    죽음 = float((np.abs(양) < 1e-12).mean())
    프레임.append([go.Heatmap(z=복원[조각][::-1], colorscale="gray",
                             zmin=0, zmax=1, showscale=False)])
    이름표.append(f"{q:g}")
배치 = dict(title=dict(text="같은 조각을 품질만 바꿔 가며 (8x8 격자를 보라)"),
           xaxis=dict(visible=False, scaleanchor="y"),
           yaxis=dict(visible=False), height=520,
           margin=dict(l=40, r=40, t=60, b=40))
slider_figure(프레임, 이름표, 배치, prefix="q = ", initial=0)

def 필요2(X, 목표=0.99):
    e = np.sort(np.asarray(X).ravel()**2)[::-1]
    return int(np.searchsorted(np.cumsum(e)/e.sum(), 목표) + 1)

Q, _ = np.linalg.qr(rng.normal(size=(m, m)))     # 아무 직교행렬
가우스 = rng.standard_normal((m, n))             # 정리가 말하는 그 백색잡음
print("정규분포 백색잡음을 여러 정규직교 기저로 옮겨 본다")
print(f"{'기저':>14}{'표준편차':>11}{'99% 에 필요한 계수':>22}")
for 이름, 계수 in (("픽셀 (원본)", 가우스),
                  ("DCT", dctn(가우스, norm="ortho")),
                  ("하르", H @ 가우스 @ H.T),
                  ("무작위 직교", Q.T @ 가우스 @ Q)):
    k = 필요2(계수)
    print(f"{이름:>14}{계수.std():>11.4f}{f'{k:,} ({k/가우스.size:.1%})':>22}")
print()
print("-> 전부 73% 대로 같다. 직교변환은 백색잡음의 분포를 바꾸지 못한다.")

잡음 = rng.random((m, n))                        # 균등분포
가운데 = 잡음 - 잡음.mean()
print(f"{'기저':>14}{'왜도':>9}{'첨도':>9}{'99% 에 필요한 계수':>22}")
for 이름, 계수 in (("픽셀 (균등)", 가운데),
                  ("DCT", dctn(가운데, norm="ortho")),
                  ("하르", H @ 가운데 @ H.T)):
    v = 계수.ravel()
    왜도 = float(((v - v.mean())**3).mean() / v.std()**3)
    첨도 = float(((v - v.mean())**4).mean() / v.std()**4)
    k = 필요2(계수)
    print(f"{이름:>14}{왜도:>9.3f}{첨도:>9.3f}{f'{k:,} ({k/잡음.size:.1%})':>22}")
print()
print("-> 균등분포의 첨도는 1.8, 정규분포는 3.0 이다.")
print("   변환을 거치면 많은 값을 더하게 되므로 중심극한정리로 정규분포에 가까워진다.")
print("   에너지를 몰아주지는 못하지만 분포 모양은 바꾼 셈이다.")

print("사진과 잡음을 견준다")
print(f"{'':>12}{'이웃 픽셀 상관':>14}{'99% 에 필요한 DCT 계수':>26}")
for 이름, X in (("시험 사진", 사진), ("무작위 잡음", 잡음)):
    a, b = X[:, :-1].ravel(), X[:, 1:].ravel()
    k = 필요2(dctn(X, norm="ortho"))
    print(f"{이름:>12}{np.corrcoef(a, b)[0,1]:>14.4f}"
          f"{f'{k:,} ({k/X.size:.1%})':>26}")
print()
print("압축이 가능한 것은 수학 덕분이 아니라 세계가 매끄럽기 때문이다.")

그림 = go.Figure(layout=dict(
    title=dict(text="DCT 계수를 큰 것부터 늘어놓으면"),
    xaxis=dict(title=dict(text="몇 번째"), type="log"),
    yaxis=dict(title=dict(text="크기"), type="log"),
    height=420, margin=dict(l=70, r=30, t=60, b=50)))
for 이름, X, 색 in (("시험 사진", 사진, COLORS["output"]),
                   ("무작위 잡음", 잡음, "#999999")):
    cc = np.sort(np.abs(dctn(X, norm="ortho")).ravel())[::-1]
    그림.add_scatter(x=np.arange(1, cc.size + 1), y=cc, mode="lines",
                    name=이름, line=dict(color=색, width=2.4))
그림