DOI: 10.5281/zenodo.17685434
title: vCF-1 Demo: Linear Curvature-Change Denoising for 2D Astronomical Maps (Medley 2025)

description: vCF-1 (Variable Curvature Filter – Version 1) is a linear, prior-minimal denoising method developed by Michael Medley as part of the Echo Lattice Theory Family v2.

This record provides a safe, simplified demonstration of how the vCF-1 curvature-change operator can be applied to 2D scientific images such as comet frames, sky maps, and astronomical data products.

 

This demo does not include the full proprietary vCF-1 implementation.

Instead, it shows a transparent, reproducible subset of the method using:

 



1D first-order and second-order finite-difference operators

A 2D curvature-change operator constructed via Kronecker products

A linear solve of the form(I + λE) x = ywhere E is the curvature-change operator and λ controls smoothing strength

A synthetic test image (Gaussian source + noise) to illustrate performance


 

 

 

What This Demo Shows

 

 



How vCF-1 removes high-frequency noise while maintaining large-scale structures

How the method applies to 2D arrays (e.g., FITS-loaded maps)

How the curvature-change framework differs from traditional wavelet or Gaussian filters

A fully reproducible workflow using only NumPy, SciPy, and Matplotlib


 

 

 

What This Demo Does Not Reveal

 

 



No full Echo Lattice engine

No production-level 2D operator

No HEALPix or domain-specific CMB code

No proprietary optimizations


 
description: #!/usr/bin/env python3# ============================================================# vCF-1 on 2D Astronomical Maps (DEMO VERSION)# © 2025 Michael Medley | Echo Lattice Theory Family v2# ------------------------------------------------------------# NOTE:# This is a *simplified* 2D demonstration of the vCF-1 idea# for images/maps (e.g., sky maps, comet frames).# It does NOT expose the full proprietary vCF-1 implementation.# ============================================================

import numpy as npimport matplotlib.pyplot as pltfrom scipy.sparse import diags, eye, kronfrom scipy.sparse.linalg import cg

# ------------------------------------------------------------# 1D DIFFERENCE OPERATORS (DEMO)# ------------------------------------------------------------

def build_D1_1d(N: int):    """1D first-difference operator (N-1 x N)."""    data = np.array([-np.ones(N - 1), np.ones(N - 1)])    offsets = np.array([0, 1])    return diags(data, offsets, shape=(N - 1, N))

def build_D2_1d(N: int):    """1D second-difference operator (N-2 x N)."""    data = np.array([np.ones(N - 2), -2 * np.ones(N - 2), np.ones(N - 2)])    offsets = np.array([0, 1, 2])    return diags(data, offsets, shape=(N - 2, N))

# ------------------------------------------------------------# 2D CURVATURE-CHANGE OPERATOR (DEMO)# ------------------------------------------------------------

def build_E_2d(nx: int, ny: int, alpha: float = 0.52):    """    Build a simple 2D curvature-change operator E for an (ny x nx) image.    We combine x- and y-direction operators using Kronecker products.    DEMO ONLY – does not represent the full Echo Lattice implementation.    """    # 1D operators in x and y    D1x = build_D1_1d(nx)    D2x = build_D2_1d(nx)    D1y = build_D1_1d(ny)    D2y = build_D2_1d(ny)

    Ix = eye(nx, format="csc")    Iy = eye(ny, format="csc")

    # 2D first differences: along x and y    J1x = kron(Iy, D1x)    # horizontal differences    J1y = kron(D1y, Ix)    # vertical differences

    # 2D second differences: along x and y    J2x = kron(Iy, D2x)    J2y = kron(D2y, Ix)

    # Stack all "curvature" rows (demo approximation)    # and build E = α J1ᵀJ1 + (1-α) J2ᵀJ2    # (Here we just sum contributions from x and y directions.)    E1 = J1x.T @ J1x + J1y.T @ J1y    E2 = J2x.T @ J2x + J2y.T @ J2y

    E = alpha * E1 + (1.0 - alpha) * E2    return E.tocsc()

def build_A_2d(nx: int, ny: int, lam: float, alpha: float = 0.52,               jitter: float = 1e-12):    """    System matrix A_2D = I + λE_2D with small diagonal jitter.    DEMO ONLY.    """    N = nx * ny    E = build_E_2d(nx, ny, alpha=alpha)    I = eye(N, format="csc")    A = I + lam * E

    # Add jitter for numerical stability    A = A.tolil()    A.setdiag(A.diagonal() + jitter)    return A.tocsc()

# ------------------------------------------------------------# 2D DENOISING (DEMO)# ------------------------------------------------------------

def vcf1_denoise_2d(img: np.ndarray,                    alpha: float = 0.52,                    lam: float = 9.8e-3):    """    Apply a vCF-1-style curvature-change filter to a 2D image (DEMO).    img is a 2D numpy array (ny x nx).

    This function flattens the image, applies (I + λE)^{-1},    then reshapes back to (ny x nx).    """    img = np.asarray(img, dtype=float)    ny, nx = img.shape    N = nx * ny

    A = build_A_2d(nx, ny, lam=lam, alpha=alpha)    y = img.reshape(N)

    x, info = cg(A, y, maxiter=2000, tol=1e-6)    if info != 0:        print(f"CG did not fully converge (info={info}), returning input.")        return img

    return x.reshape((ny, nx))

# ------------------------------------------------------------# DEMO DRIVER: SYNTHETIC IMAGE# ------------------------------------------------------------

def demo_synthetic():    """    Demonstration on a synthetic 'astronomical' image:    a smooth Gaussian blob + noise.    """    ny, nx = 128, 128    y_grid, x_grid = np.mgrid[0:ny, 0:nx]    cx, cy = nx / 2, ny / 2    r2 = (x_grid - cx)**2 + (y_grid - cy)**2

    # Smooth "source" + noise    img_clean_true = np.exp(-r2 / (2 * (0.15 * nx)**2))    rng = np.random.default_rng(42)    noise = 0.2 * rng.standard_normal(size=(ny, nx))    img_noisy = img_clean_true + noise

    img_vcf1 = vcf1_denoise_2d(img_noisy, alpha=0.52, lam=9.8e-3)

    # Plot before/after    vmin = np.min(img_noisy)    vmax = np.max(img_noisy)

    fig, axs = plt.subplots(1, 3, figsize=(12, 4))    im0 = axs[0].imshow(img_noisy, origin="lower", vmin=vmin, vmax=vmax)    axs[0].set_title("Noisy input (demo)")    plt.colorbar(im0, ax=axs[0], fraction=0.046, pad=0.04)

    im1 = axs[1].imshow(img_vcf1, origin="lower", vmin=vmin, vmax=vmax)    axs[1].set_title("vCF-1 cleaned (demo)")    plt.colorbar(im1, ax=axs[1], fraction=0.046, pad=0.04)

    im2 = axs[2].imshow(img_clean_true, origin="lower", vmin=vmin, vmax=vmax)    axs[2].set_title("True underlying signal")    plt.colorbar(im2, ax=axs[2], fraction=0.046, pad=0.04)

    for ax in axs:        ax.set_xticks([])        ax.set_yticks([])

    plt.tight_layout()    plt.show()

if __name__ == "__main__":    # Run the synthetic demo.    # For real data, replace img_noisy with any 2D map (e.g. FITS-loaded map).    demo_synthetic()