diff --git a/Source/cielim-python/cielim/image_comparison_toolkit.py b/Source/cielim-python/cielim/image_comparison_toolkit.py
new file mode 100644
index 00000000..1e4a8883
--- /dev/null
+++ b/Source/cielim-python/cielim/image_comparison_toolkit.py
@@ -0,0 +1,226 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from PIL import Image
+
+THRESHOLD = 10
+
+
+def load_grayscale(path):
+    return np.array(Image.open(path).convert("L"))
+
+
+def compute_disk_stats(img, threshold=THRESHOLD):
+    pixels = img[img > threshold]
+    return {
+        "pixels": pixels,
+        "mean": np.mean(pixels),
+        "std": np.std(pixels),
+    }
+
+
+def plot_histogram(img1, img2, title1="img1", title2="img2", bins=256):
+    s1 = compute_disk_stats(img1)
+    s2 = compute_disk_stats(img2)
+
+    bin_edges = np.linspace(0, 255, bins + 1)
+    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+    counts1, _ = np.histogram(img1.flatten(), bins=bins, range=(0, 255))
+    counts2, _ = np.histogram(img2.flatten(), bins=bins, range=(0, 255))
+
+    fig, ax = plt.subplots(figsize=(12, 6))
+    fig.suptitle("Pixel Intensity Distribution & Disk Statistics", fontsize=14, fontweight="bold")
+
+    ax.fill_between(bin_centers, counts1, alpha=0.4, color="red", label=title1)
+    ax.fill_between(bin_centers, counts2, alpha=0.4, color="blue", label=title2)
+    ax.step(bin_centers, counts1, color="darkred", alpha=0.9, linewidth=1.0)
+    ax.step(bin_centers, counts2, color="darkblue", alpha=0.9, linewidth=1.0)
+
+    ax.axvline(s1["mean"], color="darkred", linewidth=2, linestyle="--", label=f"{title1} mean: {s1['mean']:.1f}")
+    ax.axvline(s2["mean"], color="darkblue", linewidth=2, linestyle="--", label=f"{title2} mean: {s2['mean']:.1f}")
+    ax.axvspan(
+        s1["mean"] - s1["std"], s1["mean"] + s1["std"], alpha=0.1, color="red", label=f"{title1} ±1σ: {s1['std']:.1f}"
+    )
+    ax.axvspan(
+        s2["mean"] - s2["std"], s2["mean"] + s2["std"], alpha=0.1, color="blue", label=f"{title2} ±1σ: {s2['std']:.1f}"
+    )
+
+    ax.set_yscale("log")
+    ax.set_xlabel("Pixel Intensity (0–255)", fontsize=12)
+    ax.set_ylabel("Number of Pixels (log scale)", fontsize=12)
+    ax.set_xlim(0, 255)
+    ax.legend(fontsize=10)
+
+    plt.tight_layout()
+    return fig
+
+
+def plot_difference(img1, img2, title1="img1", title2="img2", bins=256):
+    s1 = compute_disk_stats(img1)
+    s2 = compute_disk_stats(img2)
+
+    bin_edges = np.linspace(0, 255, bins + 1)
+    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+    counts1, _ = np.histogram(s1["pixels"], bins=bins, range=(0, 255))
+    counts2, _ = np.histogram(s2["pixels"], bins=bins, range=(0, 255))
+    diff = counts1.astype(int) - counts2.astype(int)
+    pos_mask = diff >= 0
+    neg_mask = diff < 0
+
+    fig, ax = plt.subplots(figsize=(12, 6))
+    fig.suptitle(f"Histogram Difference — {title1} vs {title2}", fontsize=14, fontweight="bold")
+
+    ax.bar(bin_centers[pos_mask], diff[pos_mask], width=1.0, color="red", alpha=0.8, label=f"{title1} > {title2}")
+    ax.bar(bin_centers[neg_mask], diff[neg_mask], width=1.0, color="blue", alpha=0.8, label=f"{title2} > {title1}")
+    ax.axhline(0, color="black", linewidth=0.8)
+    ax.set_xlabel("Pixel Intensity (0–255)", fontsize=12)
+    ax.set_ylabel(f"Count Difference ({title1} − {title2})", fontsize=12)
+    ax.set_xlim(0, 255)
+    ax.legend(fontsize=10)
+
+    plt.tight_layout()
+    return fig
+
+
+def match_shapes(img1, img2):
+    if img1.shape != img2.shape:
+        img2 = np.array(Image.fromarray(img2).resize((img1.shape[1], img1.shape[0]), Image.BILINEAR))
+    return img1, img2
+
+
+def plot_diff_heatmap(img1, img2, title1="img1", title2="img2"):
+    img1, img2 = match_shapes(img1, img2)
+    mask = (img1 > THRESHOLD) | (img2 > THRESHOLD)
+    diff = np.zeros_like(img1, dtype=float)
+    diff[mask] = img1[mask].astype(float) - img2[mask].astype(float)
+
+    fig, ax = plt.subplots(figsize=(8, 8))
+    fig.suptitle(f"Pixel Difference Heatmap — {title1} minus {title2}", fontsize=13, fontweight="bold")
+
+    im = ax.imshow(diff, cmap="RdBu", vmin=-128, vmax=128)
+    cbar = plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
+    cbar.set_label(f"Intensity Difference ({title1} − {title2})", fontsize=10)
+
+    ax.set_title(f"Red = {title1} brighter  |  Blue = {title2} brighter  |  White = no difference", fontsize=9)
+    ax.axis("off")
+
+    plt.tight_layout()
+    return fig
+
+
+def generate_plots(img1, img2, output_dir, title1="img1", title2="img2", align=False, reference=1):
+    """
+    reference=1 means img1 is the reference for alignment (default).
+    reference=2 means img2 is the reference.
+    """
+    from pathlib import Path
+    import numpy as np
+
+    output_dir = Path(output_dir)
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    if not isinstance(img1, np.ndarray):
+        img1 = load_grayscale(img1)
+    if not isinstance(img2, np.ndarray):
+        img2 = load_grayscale(img2)
+
+    fig1 = plot_histogram(img1, img2, title1=title1, title2=title2)
+    fig1.savefig(output_dir / "histogram.png", dpi=150, bbox_inches="tight")
+    plt.close(fig1)
+
+    fig2 = plot_difference(img1, img2, title1=title1, title2=title2)
+    fig2.savefig(output_dir / "histogram_difference.png", dpi=150, bbox_inches="tight")
+    plt.close(fig2)
+
+    fig3 = plot_diff_heatmap(img1, img2, title1=title1, title2=title2)
+    fig3.savefig(output_dir / "heatmap_diff.png", dpi=150, bbox_inches="tight")
+    plt.close(fig3)
+
+    if align:
+        if reference == 1:
+            fig4, shift_y, shift_x = plot_alignment(img1, img2, title1=title1, title2=title2)
+        else:
+            fig4, shift_y, shift_x = plot_alignment(img2, img1, title1=title2, title2=title1)
+        fig4.savefig(output_dir / "alignment.png", dpi=150, bbox_inches="tight")
+        plt.close(fig4)
+
+
+def cross_correlate_fft(img1, img2):
+    i1 = img1.astype(float) - img1.mean()
+    i2 = img2.astype(float) - img2.mean()
+    corr = np.fft.ifftshift(np.real(np.fft.ifft2(np.fft.fft2(i1) * np.conj(np.fft.fft2(i2)))))
+    peak = np.unravel_index(np.argmax(corr), corr.shape)
+    shift_y = peak[0] - corr.shape[0] // 2
+    shift_x = peak[1] - corr.shape[1] // 2
+    return shift_y, shift_x, corr
+
+
+def apply_shift(img, shift_y, shift_x):
+    return np.roll(np.roll(img, shift_y, axis=0), shift_x, axis=1)
+
+
+def plot_alignment(img1, img2, title1="img1", title2="img2"):
+    img1, img2 = match_shapes(img1, img2)
+
+    shift_y, shift_x, corr = cross_correlate_fft(img1, img2)
+    img2_aligned = apply_shift(img2, shift_y, shift_x)
+
+    mask_b = (img1 > THRESHOLD) | (img2 > THRESHOLD)
+    diff_before = np.zeros_like(img1, dtype=float)
+    diff_before[mask_b] = img1[mask_b].astype(float) - img2[mask_b].astype(float)
+
+    mask_a = (img1 > THRESHOLD) | (img2_aligned > THRESHOLD)
+    diff_after = np.zeros_like(img1, dtype=float)
+    diff_after[mask_a] = img1[mask_a].astype(float) - img2_aligned[mask_a].astype(float)
+
+    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
+    fig.suptitle(
+        f"Cross-Correlation Alignment\n"
+        f"Reference: {title1}  |  Aligned: {title2}  |  "
+        f"Offset — i (x): {shift_x}px,  j (y): {shift_y}px",
+        fontsize=13,
+        fontweight="bold",
+    )
+
+    axes[0, 0].imshow(img1, cmap="gray", vmin=0, vmax=255)
+    axes[0, 0].set_title(f"{title1}\n(reference)", fontsize=11)
+    axes[0, 0].axis("off")
+
+    axes[0, 1].imshow(img2, cmap="gray", vmin=0, vmax=255)
+    axes[0, 1].set_title(f"{title2} — before alignment", fontsize=11)
+    axes[0, 1].axis("off")
+
+    axes[0, 2].imshow(img2_aligned, cmap="gray", vmin=0, vmax=255)
+    axes[0, 2].set_title(f"{title2} — after alignment\n(i={shift_x}px, j={shift_y}px)", fontsize=11)
+    axes[0, 2].axis("off")
+
+    ax_vec = axes[1, 0]
+    lim = max(abs(shift_x), abs(shift_y), 10) * 1.5
+    ax_vec.set_xlim(-lim, lim)
+    ax_vec.set_ylim(-lim, lim)
+    ax_vec.axhline(0, color="gray", linewidth=0.8, linestyle="--")
+    ax_vec.axvline(0, color="gray", linewidth=0.8, linestyle="--")
+    ax_vec.annotate("", xy=(shift_x, 0), xytext=(0, 0), arrowprops=dict(arrowstyle="->", color="blue", lw=2))
+    ax_vec.annotate("", xy=(0, -shift_y), xytext=(0, 0), arrowprops=dict(arrowstyle="->", color="red", lw=2))
+    ax_vec.annotate("", xy=(shift_x, -shift_y), xytext=(0, 0), arrowprops=dict(arrowstyle="->", color="green", lw=2.5))
+    ax_vec.plot(shift_x, -shift_y, "go", markersize=8)
+    ax_vec.text(shift_x + lim * 0.05, lim * 0.05, f"i={shift_x}px", color="blue", fontsize=10)
+    ax_vec.text(lim * 0.05, -shift_y + lim * 0.05, f"j={shift_y}px", color="red", fontsize=10)
+    ax_vec.text(shift_x + lim * 0.05, -shift_y + lim * 0.05, f"({shift_x}, {shift_y})", color="green", fontsize=9)
+    ax_vec.set_title(f"Offset vector\n{title2} → {title1}", fontsize=11)
+    ax_vec.set_xlabel("i (x shift) px")
+    ax_vec.set_ylabel("j (y shift) px")
+    ax_vec.set_aspect("equal")
+    ax_vec.grid(True, alpha=0.3)
+
+    im1 = axes[1, 1].imshow(diff_before, cmap="RdBu", vmin=-128, vmax=128)
+    plt.colorbar(im1, ax=axes[1, 1], fraction=0.046, pad=0.04)
+    axes[1, 1].set_title(f"Heatmap — BEFORE alignment\n({title1} − {title2})", fontsize=11)
+    axes[1, 1].axis("off")
+
+    im2 = axes[1, 2].imshow(diff_after, cmap="RdBu", vmin=-128, vmax=128)
+    plt.colorbar(im2, ax=axes[1, 2], fraction=0.046, pad=0.04)
+    axes[1, 2].set_title(f"Heatmap — AFTER alignment\n({title1} − {title2} aligned)", fontsize=11)
+    axes[1, 2].axis("off")
+
+    plt.tight_layout()
+    return fig, shift_y, shift_x
diff --git a/Source/cielim-python/cielim/tests/test_image_comparison_toolkit.py b/Source/cielim-python/cielim/tests/test_image_comparison_toolkit.py
new file mode 100644
index 00000000..6864586c
--- /dev/null
+++ b/Source/cielim-python/cielim/tests/test_image_comparison_toolkit.py
@@ -0,0 +1,53 @@
+import numpy as np
+import pytest
+from pathlib import Path
+import sys
+
+HERE = Path(__file__).resolve().parent
+CIELIM_ROOT = HERE.parent
+
+sys.path.insert(0, str(CIELIM_ROOT))
+
+from image_comparison_toolkit import compute_disk_stats, cross_correlate_fft, THRESHOLD
+
+
+@pytest.fixture
+def circle_image():
+
+    size = 100
+    img = np.zeros((size, size), dtype=np.uint8)
+    cy, cx = size // 2, size // 2
+    y, x = np.ogrid[:size, :size]
+    mask = (x - cx) ** 2 + (y - cy) ** 2 <= 10**2
+    img[mask] = 255
+    return img
+
+
+def test_histogram_diff_identical(circle_image):
+
+    s1 = compute_disk_stats(circle_image)
+    s2 = compute_disk_stats(circle_image)
+
+    np.testing.assert_allclose(s1["mean"], s2["mean"], rtol=0, atol=0, err_msg="Means differ for identical images")
+    np.testing.assert_allclose(
+        s1["std"], s2["std"], rtol=0, atol=0, err_msg="Standard deviations differ for identical images"
+    )
+
+    counts1, _ = np.histogram(s1["pixels"], bins=256, range=(0, 255))
+    counts2, _ = np.histogram(s2["pixels"], bins=256, range=(0, 255))
+    diff = counts1.astype(int) - counts2.astype(int)
+
+    np.testing.assert_allclose(
+        diff, np.zeros_like(diff), rtol=0, atol=0, err_msg="Expected zero histogram difference for identical images"
+    )
+
+
+@pytest.mark.parametrize("correlate_fn", [cross_correlate_fft])
+def test_cross_correlation_function(circle_image, correlate_fn):
+
+    shifted = np.roll(np.roll(circle_image, shift=20, axis=1), shift=15, axis=0)
+    shift_y, shift_x, _ = correlate_fn(circle_image, shifted)
+
+    np.testing.assert_allclose(
+        [-20, -15], [shift_x, shift_y], rtol=0, atol=0, err_msg="Cross-correlation did not recover the correct shift"
+    )
diff --git a/Source/cielim-python/examples/cielim_comparison.py b/Source/cielim-python/examples/cielim_comparison.py
new file mode 100644
index 00000000..deb0a324
--- /dev/null
+++ b/Source/cielim-python/examples/cielim_comparison.py
@@ -0,0 +1,96 @@
+from pathlib import Path
+import sys
+import cv2
+import numpy as np
+from PIL import Image
+
+HERE = Path(__file__).resolve().parent
+CIELIM_ROOT = HERE.parent
+
+sys.path.insert(0, str(CIELIM_ROOT))
+
+import context
+from image_comparison_toolkit import generate_plots
+
+BASE_DIR = CIELIM_ROOT / "support-data" / "giant-vesta"
+
+PAD_SMALL = 1.4
+PAD_LARGE = 0.20
+SMALL_THRESHOLD = 50
+DISPLAY_SIZE = 300
+ALIGN = True  # flag to enable cross-correlation alignment
+REFERENCE = 1  # 1 = img1 (cielim) is reference, 2 = img2 (giant) is reference
+
+SESSIONS = [1, 2, 3]
+
+
+def load_grayscale(path):
+    return np.array(Image.open(path).convert("L"))
+
+
+def get_cob_and_bbox(img):
+    thresh_val = max(img.max() // 4, 10)
+    _, thresh = cv2.threshold(img, thresh_val, 255, cv2.THRESH_BINARY)
+    contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
+    if not contours:
+        return img.shape[1] // 2, img.shape[0] // 2, 0, 0
+    largest = max(contours, key=cv2.contourArea)
+    bx, by, bw, bh = cv2.boundingRect(largest)
+    m = cv2.moments(thresh)
+    cx = int(m["m10"] / m["m00"]) if m["m00"] > 0 else bx + bw // 2
+    cy = int(m["m01"] / m["m00"]) if m["m00"] > 0 else by + bh // 2
+    return cx, cy, bw, bh
+
+
+def detect_pad(img):
+    _, _, bw, bh = get_cob_and_bbox(img)
+    obj_size = max(bw, bh)
+    pad = PAD_SMALL if obj_size < SMALL_THRESHOLD else PAD_LARGE
+    return pad, obj_size
+
+
+def crop_roi(img, pad_frac):
+    cx, cy, bw, bh = get_cob_and_bbox(img)
+    pad = int(max(bw, bh) * pad_frac)
+    half = max(bw, bh) // 2 + pad
+    x1 = max(cx - half, 0)
+    x2 = min(cx + half, img.shape[1])
+    y1 = max(cy - half, 0)
+    y2 = min(cy + half, img.shape[0])
+    crop = img[y1:y2, x1:x2]
+    if max(crop.shape) < DISPLAY_SIZE:
+        crop = np.array(Image.fromarray(crop).resize((DISPLAY_SIZE, DISPLAY_SIZE), Image.NEAREST))
+    return crop
+
+
+if __name__ == "__main__":
+    for session in SESSIONS:
+        cielim_img = BASE_DIR / f"cielim_{session}.png"
+        giant_img = BASE_DIR / f"giant_{session}.png"
+        out_dir = BASE_DIR / "comparison_plots" / f"session_{session}"
+
+        if not cielim_img.exists():
+            raise FileNotFoundError(f"Missing: {cielim_img}")
+        if not giant_img.exists():
+            raise FileNotFoundError(f"Missing: {giant_img}")
+
+        img_c = load_grayscale(cielim_img)
+        img_g = load_grayscale(giant_img)
+
+        pad_c, _ = detect_pad(img_c)
+        pad_g, _ = detect_pad(img_g)
+
+        crop_c = crop_roi(img_c, pad_c)
+        crop_g = crop_roi(img_g, pad_g)
+
+        generate_plots(
+            crop_c,
+            crop_g,
+            title1=cielim_img.stem,
+            title2=giant_img.stem,
+            output_dir=out_dir,
+            align=ALIGN,
+            reference=REFERENCE,
+        )
+
+        print(f"Session {session} done → {out_dir}")