"""Literal Python conversion of fig10_15.m."""

from __future__ import annotations

import time
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image


BASE_DIR = Path(__file__).resolve().parent


def fig10_15() -> None:
    # image separation: includes figure 10.16

    start_time = time.perf_counter()

    # read in original images
    x1 = Image.open(BASE_DIR / "flowers.jpg").convert("L")
    I1 = np.asarray(x1, dtype=float)
    size1 = I1.shape
    print("size1 =", size1)

    x2 = Image.open(BASE_DIR / "fox2.jpg").convert("L")
    I2B = np.asarray(x2, dtype=float)
    I2 = I2B[0:720, :]
    size2 = I2.shape
    print("size2 =", size2)

    fig1 = plt.figure(1, figsize=(8.0, 4.5))
    fig1.clf()
    ax1 = fig1.add_subplot(111)
    ax1.imshow(I1, cmap="gray", vmin=np.min(I1), vmax=np.max(I1))
    # ax1.imagesc(I1)
    ax1.set_aspect("equal")
    ax1.axis("off")
    fig1.tight_layout(pad=0)

    fig2 = plt.figure(2, figsize=(8.0, 4.5))
    fig2.clf()
    ax2 = fig2.add_subplot(111)
    ax2.imshow(I2, cmap="gray", vmin=np.min(I2), vmax=np.max(I2))
    # ax2.imagesc(I2)
    ax2.set_aspect("equal")
    ax2.axis("off")
    fig2.tight_layout(pad=0)

    MM = size1[0]
    NN = size1[1]
    n = MM * NN

    # mixing matrix
    M11 = 0.75
    M12 = 0.25
    M21 = 0.9
    M22 = 0.1
    M = np.array([[M11, M12], [M21, M22]])

    # calculate and show the mixed images
    III1 = M11 * I1 + M12 * I2
    III2 = M21 * I1 + M22 * I2

    fig3 = plt.figure(3, figsize=(8.0, 4.5))
    fig3.clf()
    ax3 = fig3.add_subplot(111)
    ax3.imshow(III1, cmap="gray", vmin=np.min(III1), vmax=np.max(III1))
    # ax3.imagesc(III1)
    ax3.set_aspect("equal")
    ax3.axis("off")
    fig3.tight_layout(pad=0)

    fig4 = plt.figure(4, figsize=(8.0, 4.5))
    fig4.clf()
    ax4 = fig4.add_subplot(111)
    ax4.imshow(III2, cmap="gray", vmin=np.min(III2), vmax=np.max(III2))
    # ax4.imagesc(III2)
    ax4.set_aspect("equal")
    ax4.axis("off")
    fig4.tight_layout(pad=0)

    # center image data
    xM1 = np.sum(III1) / n
    xM2 = np.sum(III2) / n
    X1 = III1 - xM1
    X2 = III2 - xM2

    # compute product matrix
    X11 = np.sum(X1 * X1) / n
    X12 = np.sum(X1 * X2) / n
    X22 = np.sum(X2 * X2) / n

    XX = np.array([[X11, X12], [X12, X22]])
    print("XX =\n", XX)

    # find factorization of inner product matrix
    eigenvalues, M1 = np.linalg.eig(XX)
    D1 = np.diag(eigenvalues)

    if D1[0, 0] >= D1[1, 1]:
        d1 = D1[0, 0]
        d2 = D1[1, 1]
        Q1 = M1[:, 0] / np.linalg.norm(M1[:, 0], 2)
        Q2 = M1[:, 1] / np.linalg.norm(M1[:, 1], 2)
    else:
        d2 = D1[0, 0]
        d1 = D1[1, 1]
        Q2 = M1[:, 0] / np.linalg.norm(M1[:, 0], 2)
        Q1 = M1[:, 1] / np.linalg.norm(M1[:, 1], 2)

    Q = np.column_stack((Q1, Q2))
    D5 = np.zeros((2, 2))
    D5[0, 0] = 1 / np.sqrt(d1)
    D5[1, 1] = 1 / np.sqrt(d2)
    QD = D5 @ Q.T

    # find y_i's
    Y1 = QD[0, 0] * X1 + QD[0, 1] * X2
    Y2 = QD[1, 0] * X1 + QD[1, 1] * X2

    # compute kurtosis as function of angle
    ntheta = 360 * 4
    print("ntheta =", ntheta)
    # ntheta=360*16
    theta = np.linspace(0, 2 * np.pi, ntheta)
    # The following vectorized fourth-moment calculation is algebraically
    # identical to the original loop, but avoids constructing two full
    # 720-by-1280 images for every value of theta.
    m40 = np.sum(Y1**4) / n
    m31 = np.sum((Y1**3) * Y2) / n
    m22 = np.sum((Y1**2) * (Y2**2)) / n
    m13 = np.sum(Y1 * (Y2**3)) / n
    m04 = np.sum(Y2**4) / n
    c = np.cos(theta)
    s = np.sin(theta)
    KKs1 = (
        c**4 * m40
        - 4 * c**3 * s * m31
        + 6 * c**2 * s**2 * m22
        - 4 * c * s**3 * m13
        + s**4 * m04
    )
    KKs2 = (
        s**4 * m40
        + 4 * s**3 * c * m31
        + 6 * s**2 * c**2 * m22
        + 4 * s * c**3 * m13
        + c**4 * m04
    )
    K = np.abs(KKs1 - 3) + np.abs(KKs2 - 3)
    # K=(KKs1-3)**2+(KKs2-3)**2

    # locate max points and then pick one for separation
    thetaMM_list: list[float] = []
    for ia in range(1, ntheta - 1):
        if K[ia - 1] < K[ia] and K[ia + 1] < K[ia]:
            thetaMM_list.append(theta[ia])
    thetaMM = np.asarray(thetaMM_list)
    print("thetaMM =", thetaMM)

    # theta0=thetaMM[1]
    theta0 = thetaMM[3]
    print("theta0 =", theta0)

    fig5 = plt.figure(5, figsize=(7.2, 3.4))
    fig5.clf()

    co = np.array(
        [
            [0, 0, 1],
            [0, 0.5, 0],
            [1, 0, 0],
            [0, 0.75, 0.75],
            [0.75, 0, 0.75],
            [0.75, 0.75, 0],
            [0.25, 0.25, 0.25],
        ]
    )
    ax5 = fig5.add_subplot(111)
    ax5.set_prop_cycle(color=co.tolist())

    ax5.plot(theta, K, linewidth=1.4)

    # ax5.axis([0, 2*np.pi, 1.5, 3])
    ax5.set_xlabel(r"$\theta$-axis", fontsize=14, fontweight="bold")
    ax5.set_ylabel("K -axis", fontsize=14, fontweight="bold")
    ax5.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])
    ax5.set_xticklabels(["0", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"])
    ax5.grid(True)
    ax5.set_frame_on(True)
    ax5.tick_params(labelsize=14)  # MATLAB FontSize 16 reduced by 2
    for label in ax5.get_xticklabels() + ax5.get_yticklabels():
        label.set_fontweight("bold")
    fig5.tight_layout()

    # find mixing matrix and then separate images
    V = np.array(
        [
            [np.cos(theta0), -np.sin(theta0)],
            [np.sin(theta0), np.cos(theta0)],
        ]
    )
    DD = np.zeros((2, 2))
    DD[0, 0] = np.sqrt(d1)
    DD[1, 1] = np.sqrt(d2)
    A = Q @ DD @ V.T
    print("A =\n", A)

    II1 = V[0, 0] * Y1 + V[0, 1] * Y2
    II2 = V[1, 0] * Y1 + V[1, 1] * Y2

    fig6 = plt.figure(6, figsize=(8.0, 4.5))
    fig6.clf()
    ax6 = fig6.add_subplot(111)
    ax6.imshow(II1, cmap="gray", vmin=np.min(II1), vmax=np.max(II1))
    # ax6.imagesc(II1)
    ax6.set_aspect("equal")
    ax6.axis("off")
    fig6.tight_layout(pad=0)

    fig7 = plt.figure(7, figsize=(8.0, 4.5))
    fig7.clf()
    ax7 = fig7.add_subplot(111)
    ax7.imshow(II2, cmap="gray", vmin=np.min(II2), vmax=np.max(II2))
    # ax7.imagesc(II2)
    ax7.set_aspect("equal")
    ax7.axis("off")
    fig7.tight_layout(pad=0)

    print(f"Elapsed time is {time.perf_counter() - start_time:.6f} seconds.")
    plt.show()


if __name__ == "__main__":
    fig10_15()
