{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "171e14ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%pip install colormath\n",
    "#%pip install numpy==1.22.0\n",
    "print(np.__version__)\n",
    "#%pip install numpy==1.22.0 matplotlib>=3.5\n",
    "import os\n",
    "print(\"Current Working Directory:\", os.getcwd())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93893f5a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from colormath.color_objects import LabColor, sRGBColor\n",
    "from colormath.color_diff import delta_e_cie2000\n",
    "from colormath.color_conversions import convert_color\n",
    "\n",
    "# ================================================\n",
    "# Helper Function\n",
    "# ================================================\n",
    "def lab_to_rgb(lab_color):\n",
    "    \"\"\"Convert Lab color to RGB (clipped to [0,1]).\"\"\"\n",
    "    try:\n",
    "        rgb = convert_color(lab_color, sRGBColor).get_value_tuple()\n",
    "        return np.clip(rgb, 0.0, 1.0)\n",
    "    except:\n",
    "        return (0, 0, 0)  # Black for out-of-gamut colors\n",
    "\n",
    "# ================================================\n",
    "# Grid Parameters\n",
    "# ================================================\n",
    "L = 50.0               # Fixed lightness\n",
    "a_min, a_max = -128, 128 # Chromatic a* range\n",
    "b_min, b_max = -128, 128 # Chromatic b* range\n",
    "step = 1               # Grid spacing\n",
    "neighbor_step = step / 1000  # Distance to neighbors\n",
    "save_plots = True      # Set to True to save plots as images\n",
    "\n",
    "a_values = np.arange(a_min, a_max + step, step)\n",
    "b_values = np.arange(b_min, b_max + step, step)\n",
    "aa, bb = np.meshgrid(a_values, b_values, indexing='ij')\n",
    "\n",
    "# 8-neighbor offsets (da, db)\n",
    "neighbor_offsets = [\n",
    "    (1, 0), (1, 1), (0, 1), (-1, 1),\n",
    "    (-1, 0), (-1, -1), (0, -1), (1, -1)\n",
    "]\n",
    "\n",
    "# ================================================\n",
    "# Precompute ΔE2000 distances to all neighbors\n",
    "# ================================================\n",
    "distances = np.full(aa.shape + (8,), np.nan)\n",
    "for i in range(aa.shape[0]):\n",
    "    for j in range(aa.shape[1]):\n",
    "        current_lab = (L, aa[i, j], bb[i, j])\n",
    "        for k, (da, db) in enumerate(neighbor_offsets):\n",
    "            a_neighbor = aa[i, j] + da * neighbor_step\n",
    "            b_neighbor = bb[i, j] + db * neighbor_step\n",
    "            if (a_min <= a_neighbor <= a_max) and (b_min <= b_neighbor <= b_max):\n",
    "                neighbor_lab = (L, a_neighbor, b_neighbor)\n",
    "                color1 = LabColor(lab_l=current_lab[0], lab_a=current_lab[1], lab_b=current_lab[2])\n",
    "                color2 = LabColor(lab_l=neighbor_lab[0], lab_a=neighbor_lab[1], lab_b=neighbor_lab[2])\n",
    "                distances[i, j, k] = delta_e_cie2000(color1, color2)\n",
    "\n",
    "# ================================================\n",
    "# Compute Anisotropy Measures\n",
    "# ================================================\n",
    "def tensor_anisotropy(dists, offsets, step):\n",
    "    \"\"\"Structure tensor-based anisotropy.\"\"\"\n",
    "    valid = ~np.isnan(dists)\n",
    "    valid_dists = dists[valid]\n",
    "    valid_offsets = [offsets[i] for i in range(8) if valid[i]]\n",
    "    \n",
    "    if len(valid_dists) == 0:\n",
    "        return np.nan\n",
    "    \n",
    "    T = np.zeros((2, 2))\n",
    "    for dist, (da, db) in zip(valid_dists, valid_offsets):\n",
    "        v = np.array([da * step, db * step])\n",
    "        T += dist * np.outer(v, v)\n",
    "    \n",
    "    eigvals = np.linalg.eigvalsh(T)\n",
    "    if eigvals[1] + eigvals[0] == 0:\n",
    "        return 0.0\n",
    "    return (eigvals[1] - eigvals[0]) / (eigvals[1] + eigvals[0])\n",
    "\n",
    "def fourier_anisotropy(dists, offsets):\n",
    "    \"\"\"Fourier-based anisotropy (k=2 component).\"\"\"\n",
    "    valid = ~np.isnan(dists)\n",
    "    valid_dists = dists[valid]\n",
    "    valid_offsets = [offsets[i] for i in range(8) if valid[i]]\n",
    "    \n",
    "    if len(valid_dists) < 3:\n",
    "        return np.nan\n",
    "    \n",
    "    angles = np.arctan2([db for (_, db) in valid_offsets], [da for (da, _) in valid_offsets])\n",
    "    a2 = np.mean(valid_dists * np.cos(2 * angles))\n",
    "    b2 = np.mean(valid_dists * np.sin(2 * angles))\n",
    "    return np.sqrt(a2**2 + b2**2)\n",
    "\n",
    "def statistical_anisotropy(dists):\n",
    "    \"\"\"Coefficient of variation.\"\"\"\n",
    "    valid_dists = dists[~np.isnan(dists)]\n",
    "    if len(valid_dists) < 2:\n",
    "        return np.nan\n",
    "    return np.std(valid_dists) / np.mean(valid_dists)\n",
    "\n",
    "def gradient_anisotropy(dists, offsets):\n",
    "    \"\"\"Gradient-based anisotropy (x vs. y).\"\"\"\n",
    "    right = dists[0] if not np.isnan(dists[0]) else 0\n",
    "    left = dists[4] if not np.isnan(dists[4]) else 0\n",
    "    up = dists[2] if not np.isnan(dists[2]) else 0\n",
    "    down = dists[6] if not np.isnan(dists[6]) else 0\n",
    "    \n",
    "    Gx = (right - left) / (2 * step)\n",
    "    Gy = (up - down) / (2 * step)\n",
    "    if np.abs(Gx) + np.abs(Gy) == 0:\n",
    "        return 0.0\n",
    "    return (np.abs(Gx) - np.abs(Gy)) / (np.abs(Gx) + np.abs(Gy))\n",
    "\n",
    "# Initialize anisotropy arrays\n",
    "tensor_ani = np.zeros_like(aa, dtype=float)\n",
    "fourier_ani = np.zeros_like(aa, dtype=float)\n",
    "stat_ani = np.zeros_like(aa, dtype=float)\n",
    "gradient_ani = np.zeros_like(aa, dtype=float)\n",
    "\n",
    "for i in range(aa.shape[0]):\n",
    "    for j in range(aa.shape[1]):\n",
    "        tensor_ani[i, j] = tensor_anisotropy(distances[i, j], neighbor_offsets, step)\n",
    "        fourier_ani[i, j] = fourier_anisotropy(distances[i, j], neighbor_offsets)\n",
    "        stat_ani[i, j] = statistical_anisotropy(distances[i, j])\n",
    "        gradient_ani[i, j] = gradient_anisotropy(distances[i, j], neighbor_offsets)\n",
    "\n",
    "# ================================================\n",
    "# Define rgb_grid for reference plot\n",
    "# ================================================\n",
    "rgb_grid = np.zeros((aa.shape[0], aa.shape[1], 3))\n",
    "for i in range(aa.shape[0]):\n",
    "    for j in range(aa.shape[1]):\n",
    "        lab = LabColor(L, aa[i, j], bb[i, j])\n",
    "        rgb_grid[i, j] = lab_to_rgb(lab)\n",
    "\n",
    "# ================================================\n",
    "# Plot Results\n",
    "# ================================================\n",
    "fontsize = 24  # Adjust as needed\n",
    "\n",
    "def plot_and_save(data, title, filename, cmap='viridis', vmin=None, vmax=None):\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    if data.ndim == 3:  # RGB image\n",
    "        im = plt.imshow(data, origin='lower', extent=[a_min, a_max, b_min, b_max])\n",
    "    else:  # Grayscale data\n",
    "        im = plt.imshow(data.T, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, extent=[a_min, a_max, b_min, b_max])\n",
    "    cbar = plt.colorbar(im) if data.ndim == 2 else None\n",
    "    plt.xlabel('a*', fontsize=fontsize)\n",
    "    plt.ylabel('b*', fontsize=fontsize)\n",
    "    plt.title(title, fontsize=fontsize)\n",
    "    plt.xticks(fontsize=fontsize)\n",
    "    plt.yticks(fontsize=fontsize)\n",
    "    plt.tight_layout()\n",
    "    if save_plots:\n",
    "        plt.savefig(filename, dpi=300, bbox_inches='tight')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# Reference Plot\n",
    "plot_and_save(rgb_grid, \"Reference Colors (Lab)\", \"reference_colors.png\")\n",
    "\n",
    "# Tensor-Based Anisotropy\n",
    "vmin_tensor = np.nanpercentile(tensor_ani, 1)\n",
    "vmax_tensor = np.nanpercentile(tensor_ani, 99)\n",
    "plot_and_save(tensor_ani, f\"Tensor-Based Anisotropy\\n(vmin={vmin_tensor:.2f}, vmax={vmax_tensor:.2f})\", \"tensor_anisotropy.png\", vmin=vmin_tensor, vmax=vmax_tensor)\n",
    "\n",
    "# Fourier-Based Anisotropy\n",
    "vmin_fourier = 0\n",
    "vmax_fourier = np.nanpercentile(fourier_ani, 99)\n",
    "plot_and_save(fourier_ani, f\"Fourier-Based Anisotropy\\n(vmin={vmin_fourier:.2f}, vmax={vmax_fourier:.2f})\", \"fourier_anisotropy.png\", vmin=vmin_fourier, vmax=vmax_fourier)\n",
    "\n",
    "# Statistical Anisotropy\n",
    "vmin_stat = np.nanpercentile(stat_ani, 1)\n",
    "vmax_stat = np.nanpercentile(stat_ani, 99)\n",
    "plot_and_save(stat_ani, f\"Statistical Anisotropy\\n(vmin={vmin_stat:.2f}, vmax={vmax_stat:.2f})\", \"statistical_anisotropy.png\", vmin=vmin_stat, vmax=vmax_stat)\n",
    "\n",
    "# Gradient-Based Anisotropy\n",
    "max_abs = np.nanmax(np.abs(gradient_ani))\n",
    "plot_and_save(gradient_ani, f\"Gradient-Based Anisotropy\\n(vmin={-max_abs:.2f}, vmax={max_abs:.2f})\", \"gradient_anisotropy.png\", cmap='coolwarm', vmin=-max_abs, vmax=max_abs)\n",
    "\n",
    "print(\"Experiments completed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe15e279",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from colormath.color_objects import LabColor\n",
    "from colormath.color_diff import delta_e_cie2000\n",
    "\n",
    "# ============================================================\n",
    "# Helper Functions\n",
    "# ============================================================\n",
    "\n",
    "def delta_e_wrapper(lab1, lab2):\n",
    "    # Convert Lab tuples to LabColor objects\n",
    "    color1 = LabColor(lab_l=lab1[0], lab_a=lab1[1], lab_b=lab1[2])\n",
    "    color2 = LabColor(lab_l=lab2[0], lab_a=lab2[1], lab_b=lab2[2])\n",
    "    return delta_e_cie2000(color1, color2)\n",
    "\n",
    "def is_valid_lab_color(c_lab):\n",
    "    L, a, b = c_lab\n",
    "    return -128 <= a <= 128 and -128 <= b <= 128\n",
    "\n",
    "def check_parallelogram(c_lab, delta_a=1, delta_b=1, absolute_error=False):\n",
    "    L, a, b = c_lab\n",
    "    a_plus = (L, a + delta_a, b)\n",
    "    b_plus = (L, a, b + delta_b)\n",
    "    midpoint = (L, a + delta_a / 2, b + delta_b / 2)\n",
    "    \n",
    "    d_ca = delta_e_wrapper(c_lab, a_plus)\n",
    "    d_cb = delta_e_wrapper(c_lab, b_plus)\n",
    "    d_ab = delta_e_wrapper(a_plus, b_plus)\n",
    "    d_cm = delta_e_wrapper(c_lab, midpoint)\n",
    "    \n",
    "    lhs = d_ca**2 + d_cb**2\n",
    "    rhs = 0.5 * d_ab**2 + 2 * d_cm**2\n",
    "    \n",
    "    if absolute_error:\n",
    "        error = abs(lhs - rhs)\n",
    "    else:\n",
    "        error = abs(lhs - rhs) / (lhs + rhs + 1e-10)  # Avoid division by zero\n",
    "    return error\n",
    "\n",
    "def plot_and_save(matrix, a_range, b_range, vmin, vmax, title, filename, save=False):\n",
    "    plt.figure(figsize=(10, 8))\n",
    "    im = plt.imshow(\n",
    "        matrix.T,\n",
    "        extent=[a_range.min(), a_range.max(), b_range.min(), b_range.max()],\n",
    "        origin='lower',\n",
    "        cmap='viridis',\n",
    "        aspect='auto',\n",
    "        vmin=vmin,  # Dynamic minimum\n",
    "        vmax=vmax   # Dynamic maximum\n",
    "    )\n",
    "    plt.colorbar(im, label='Error')\n",
    "    plt.xlabel('a*')\n",
    "    plt.ylabel('b*')\n",
    "    plt.title(title)\n",
    "    if save:\n",
    "        plt.savefig(filename, dpi=300, bbox_inches='tight')\n",
    "    plt.show()\n",
    "\n",
    "# ============================================================\n",
    "# Grid Parameters\n",
    "# ============================================================\n",
    "L = 50.0                # Fixed lightness\n",
    "a_min, a_max = -128, 128 # Chromatic a* range\n",
    "b_min, b_max = -128, 128 # Chromatic b* range\n",
    "step = 1                # Grid spacing\n",
    "neighbor_step = step / 2 # Distance to neighbors (default: step / 2)\n",
    "save_plots = True       # Set to True to save plots as images\n",
    "\n",
    "a_values = np.arange(a_min, a_max + step, step)\n",
    "b_values = np.arange(b_min, b_max + step, step)\n",
    "aa, bb = np.meshgrid(a_values, b_values, indexing='ij')\n",
    "\n",
    "# Variable for radius experiments\n",
    "radius_values = [0.01, 0.001]\n",
    "\n",
    "# ============================================================\n",
    "# Error Computation for Multiple Radii\n",
    "# ============================================================\n",
    "for radius in radius_values:\n",
    "    print(f\"Running experiment with radius={radius}...\")\n",
    "    error_matrix_relative = np.zeros_like(aa, dtype=np.float64)\n",
    "    error_matrix_absolute = np.zeros_like(aa, dtype=np.float64)\n",
    "    \n",
    "    for i, a in enumerate(a_values):\n",
    "        for j, b in enumerate(b_values):\n",
    "            c_lab = (L, a, b)\n",
    "            if is_valid_lab_color(c_lab):\n",
    "                try:\n",
    "                    # Compute both relative and absolute errors\n",
    "                    error_matrix_relative[i, j] = check_parallelogram(c_lab, delta_a=radius, delta_b=radius, absolute_error=False)\n",
    "                    error_matrix_absolute[i, j] = check_parallelogram(c_lab, delta_a=radius, delta_b=radius, absolute_error=True)\n",
    "                except:\n",
    "                    error_matrix_relative[i, j] = np.nan\n",
    "                    error_matrix_absolute[i, j] = np.nan\n",
    "            else:\n",
    "                error_matrix_relative[i, j] = np.nan\n",
    "                error_matrix_absolute[i, j] = np.nan\n",
    "    \n",
    "    # Replace NaNs with zeros\n",
    "    error_matrix_relative = np.nan_to_num(error_matrix_relative, nan=0.0)\n",
    "    error_matrix_absolute = np.nan_to_num(error_matrix_absolute, nan=0.0)\n",
    "    \n",
    "    # Compute dynamic range for plotting\n",
    "    vmin_rel = np.nanpercentile(error_matrix_relative, 1)\n",
    "    vmax_rel = np.nanpercentile(error_matrix_relative, 99)\n",
    "    vmin_abs = np.nanpercentile(error_matrix_absolute, 1)\n",
    "    vmax_abs = np.nanpercentile(error_matrix_absolute, 99)\n",
    "    \n",
    "    # Plot relative error\n",
    "    plot_and_save(\n",
    "        error_matrix_relative,\n",
    "        a_values,\n",
    "        b_values,\n",
    "        vmin_rel,\n",
    "        vmax_rel,\n",
    "        title=f'Relative Error (Radius={radius})',\n",
    "        filename=f'relative_error_radius_{radius}.png',\n",
    "        save=save_plots\n",
    "    )\n",
    "    \n",
    "    # Plot absolute error\n",
    "    plot_and_save(\n",
    "        error_matrix_absolute,\n",
    "        a_values,\n",
    "        b_values,\n",
    "        vmin_abs,\n",
    "        vmax_abs,\n",
    "        title=f'Absolute Error (Radius={radius})',\n",
    "        filename=f'absolute_error_radius_{radius}.png',\n",
    "        save=save_plots\n",
    "    )\n",
    "\n",
    "print(\"Experiments completed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a99e01db-c329-4711-929a-2a6b5a07ad90",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from colormath.color_objects import LabColor\n",
    "from colormath.color_diff import delta_e_cie2000\n",
    "\n",
    "# ================================================\n",
    "# Grid Parameters\n",
    "# ================================================\n",
    "L = 50.0                  # Fixed lightness\n",
    "a_min, a_max = -128, 128  # Chromatic a* range\n",
    "b_min, b_max = -128, 128  # Chromatic b* range\n",
    "step = 1                  # Grid spacing\n",
    "neighbor_step = step / 2  # Distance to neighbors\n",
    "radius_values = [0.01, 0.001]  # Radius for delta_a and delta_b\n",
    "save_plots = True         # Save generated plots\n",
    "fontsize = 24             # Font size for labels and ticks\n",
    "\n",
    "# Generate grid\n",
    "a_values = np.arange(a_min, a_max + step, step)\n",
    "b_values = np.arange(b_min, b_max + step, step)\n",
    "aa, bb = np.meshgrid(a_values, b_values, indexing='ij')\n",
    "\n",
    "# ================================================\n",
    "# Helper Functions\n",
    "# ================================================\n",
    "def delta_e_wrapper(lab1, lab2):\n",
    "    \"\"\"Compute Delta E 2000 between two Lab colors.\"\"\"\n",
    "    color1 = LabColor(lab_l=lab1[0], lab_a=lab1[1], lab_b=lab1[2])\n",
    "    color2 = LabColor(lab_l=lab2[0], lab_a=lab2[1], lab_b=lab2[2])\n",
    "    return delta_e_cie2000(color1, color2)\n",
    "\n",
    "def is_valid_lab_color(c_lab):\n",
    "    \"\"\"Check if the Lab color is within valid bounds.\"\"\"\n",
    "    _, a, b = c_lab\n",
    "    return -128 <= a <= 128 and -128 <= b <= 128\n",
    "\n",
    "def check_parallelogram(c_lab, delta_a=1, delta_b=1, neighbor_step=0.5, all_neighbors=False, absolute_error=False):\n",
    "    \"\"\"\n",
    "    Compute the error in the parallelogram law for a given CIELAB color sample.\n",
    "    \"\"\"\n",
    "    L, a, b = c_lab\n",
    "    \n",
    "    # Define all 8 neighbors in the auxiliary grid\n",
    "    neighbors = [\n",
    "        (L, a + neighbor_step, b),\n",
    "        (L, a - neighbor_step, b),\n",
    "        (L, a, b + neighbor_step),\n",
    "        (L, a, b - neighbor_step),\n",
    "        (L, a + neighbor_step, b + neighbor_step),\n",
    "        (L, a - neighbor_step, b + neighbor_step),\n",
    "        (L, a + neighbor_step, b - neighbor_step),\n",
    "        (L, a - neighbor_step, b - neighbor_step)\n",
    "    ]\n",
    "\n",
    "    # Function to compute parallelogram error for given A, B, C, D\n",
    "    def parallelogram_error(A, B, C, D):\n",
    "        d_ab = delta_e_wrapper(A, B)\n",
    "        d_bc = delta_e_wrapper(B, C)\n",
    "        d_ac = delta_e_wrapper(A, C)\n",
    "        d_bd = delta_e_wrapper(B, D)\n",
    "        \n",
    "        lhs = 2 * d_ab**2 + 2 * d_bc**2\n",
    "        rhs = d_ac**2 + d_bd**2\n",
    "        \n",
    "        if absolute_error:\n",
    "            error = abs(lhs - rhs)\n",
    "        else:\n",
    "            error = abs(lhs - rhs) / (lhs + rhs + 1e-10)  # Avoid division by zero\n",
    "        return error\n",
    "\n",
    "    if not all_neighbors:\n",
    "        # Use fixed delta_a and delta_b to define neighbors\n",
    "        B = (L, a + delta_a, b)\n",
    "        C = (L, a, b + delta_b)\n",
    "        D = (L, a + delta_a, b + delta_b)  # Form parallelogram with A, B, C\n",
    "        if is_valid_lab_color(B) and is_valid_lab_color(C) and is_valid_lab_color(D):\n",
    "            return parallelogram_error(c_lab, B, C, D)\n",
    "        else:\n",
    "            return np.nan\n",
    "    else:\n",
    "        # Iterate over all possible pairs of neighbors to find maximum error\n",
    "        max_error = 0.0\n",
    "        for i, B in enumerate(neighbors):\n",
    "            for j, C in enumerate(neighbors):\n",
    "                if i == j:  # Skip if B and C are the same\n",
    "                    continue\n",
    "                D = (L, B[1] + (C[1] - a), B[2] + (C[2] - b))  # Form D to complete parallelogram\n",
    "                if is_valid_lab_color(B) and is_valid_lab_color(C) and is_valid_lab_color(D):\n",
    "                    error = parallelogram_error(c_lab, B, C, D)\n",
    "                    max_error = max(max_error, error)\n",
    "        return max_error\n",
    "\n",
    "def plot_and_save(data, vmin, vmax, filename, title, xlabel, ylabel, fontsize=24):\n",
    "    \"\"\"Plot and optionally save the generated figure.\"\"\"\n",
    "    plt.figure(figsize=(10, 8))\n",
    "    im = plt.imshow(\n",
    "        data.T,\n",
    "        extent=[a_min, a_max, b_min, b_max],\n",
    "        origin='lower',\n",
    "        cmap='viridis',\n",
    "        aspect='auto',\n",
    "        vmin=vmin,\n",
    "        vmax=vmax\n",
    "    )\n",
    "    #plt.colorbar(im, label='Error')\n",
    "    plt.title(title, fontsize=fontsize)\n",
    "    plt.xlabel(xlabel, fontsize=fontsize)\n",
    "    plt.ylabel(ylabel, fontsize=fontsize)\n",
    "    plt.xticks(fontsize=fontsize - 4)\n",
    "    plt.yticks(fontsize=fontsize - 4)\n",
    "    if save_plots:\n",
    "        plt.savefig(filename, dpi=300, bbox_inches='tight')\n",
    "    plt.show()\n",
    "\n",
    "# ================================================\n",
    "# Main Loop: Compute and Plot Errors\n",
    "# ================================================\n",
    "for radius in radius_values:\n",
    "    delta_a = delta_b = radius\n",
    "    error_matrix_abs = np.zeros((len(a_values), len(b_values)))\n",
    "    error_matrix_rel = np.zeros((len(a_values), len(b_values)))\n",
    "    \n",
    "    # Compute errors across the grid\n",
    "    for i, a in enumerate(a_values):\n",
    "        for j, b in enumerate(b_values):\n",
    "            c_lab = (L, a, b)\n",
    "            if is_valid_lab_color(c_lab):\n",
    "                try:\n",
    "                    error = check_parallelogram(\n",
    "                        c_lab, \n",
    "                        delta_a=delta_a, \n",
    "                        delta_b=delta_b, \n",
    "                        neighbor_step=neighbor_step, \n",
    "                        all_neighbors=True,  # Test all neighbors\n",
    "                        absolute_error=False  # Use relative error\n",
    "                    )\n",
    "                    error_matrix_rel[i, j] = error\n",
    "                except:\n",
    "                    error_matrix_rel[i, j] = np.nan  # Handle edge cases\n",
    "\n",
    "                try:\n",
    "                    error = check_parallelogram(\n",
    "                        c_lab, \n",
    "                        delta_a=delta_a, \n",
    "                        delta_b=delta_b, \n",
    "                        neighbor_step=neighbor_step, \n",
    "                        all_neighbors=True,  # Test all neighbors\n",
    "                        absolute_error=True  # Use relative error\n",
    "                    )\n",
    "                    error_matrix_abs[i, j] = error\n",
    "                except:\n",
    "                    error_matrix_abs[i, j] = np.nan  # Handle edge cases\n",
    "            else:\n",
    "                error_matrix_rel[i, j] = np.nan  # Mark out-of-gamut colors as NaN\n",
    "                error_matrix_abs[i, j] = np.nan  # Mark out-of-gamut colors as NaN\n",
    "\n",
    "    # Replace NaNs with zeros\n",
    "    error_matrix_rel = np.nan_to_num(error_matrix_rel, nan=0.0)\n",
    "    error_matrix_abs = np.nan_to_num(error_matrix_abs, nan=0.0)\n",
    "\n",
    "    # Compute dynamic range for plotting\n",
    "    vmin = np.nanpercentile(error_matrix_rel, 1)\n",
    "    vmax = np.nanpercentile(error_matrix_rel, 99)\n",
    "\n",
    "    # Generate plot title and filename\n",
    "    title = f'Relative Error (spacing={radius})'\n",
    "    filename = f'parallelogram_rel_error_spacing_{radius}.png'\n",
    "    xlabel = 'a*'\n",
    "    ylabel = 'b*'\n",
    "\n",
    "    # Plot and save\n",
    "    plot_and_save(\n",
    "        data=error_matrix_rel,\n",
    "        vmin=vmin,\n",
    "        vmax=vmax,\n",
    "        filename=filename,\n",
    "        title=title,\n",
    "        xlabel=xlabel,\n",
    "        ylabel=ylabel,\n",
    "        fontsize=fontsize\n",
    "    )\n",
    "\n",
    "    # Compute dynamic range for plotting\n",
    "    vmin = np.nanpercentile(error_matrix_abs, 1)\n",
    "    vmax = np.nanpercentile(error_matrix_abs, 99)\n",
    "\n",
    "    # Generate plot title and filename\n",
    "    title = f'Absolute Error (spacing={radius})'\n",
    "    filename = f'parallelogram_abs_error_spacing_{radius}.png'\n",
    "    xlabel = 'a*'\n",
    "    ylabel = 'b*'\n",
    "\n",
    "    # Plot and save\n",
    "    plot_and_save(\n",
    "        data=error_matrix_abs,\n",
    "        vmin=vmin,\n",
    "        vmax=vmax,\n",
    "        filename=filename,\n",
    "        title=title,\n",
    "        xlabel=xlabel,\n",
    "        ylabel=ylabel,\n",
    "        fontsize=fontsize\n",
    "    )\n",
    "    \n",
    "print(\"Experiments completed\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (Color Space)",
   "language": "python",
   "name": "color_space_env"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.21"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}