impactof_Kukernel_Figures.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7492ca67-90c7-46eb-a4c8-75df747fd918",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af5d48d6-983f-4316-bc8d-fa5c224ea764",
   "metadata": {},
   "source": [
    "# Understanding Figure 2 of Gary M. Bernstein & Daniel Gruen 2014 article\n",
    "https://arxiv.org/pdf/1401.2636.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dcaeca6-a733-4d77-a750-4cc7ece88197",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.gridspec as gridspec\n",
    "params = {'legend.fontsize':14,\n",
    "          'figure.figsize': (6, 6),\n",
    "         'axes.labelsize': 16,\n",
    "         'axes.titlesize': 16,\n",
    "         'xtick.labelsize':16,\n",
    "         'ytick.labelsize': 16}\n",
    "pylab.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "181d3e47-3dd3-47d8-a043-c8802737d52e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Kernel in 1D\n",
    "from interp1DKernels import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "599e4fc8-8121-40a5-ac0b-0e527ef4b7f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Periodic kernel 1D to approx Sinc-periodic\n",
    "from periodicKernel import *  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4171e523-b2c1-44a5-b68e-b30513e0ff3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bernstein & Gruen  DFT and inverse DFT using numpy.fft code\n",
    "from BG_DFT_impl import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93a4908f-d006-4da5-a1a9-2ca30bbdd125",
   "metadata": {},
   "source": [
    "# Utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18ec0f53-2539-49c0-906d-2cce5ca5e079",
   "metadata": {},
   "outputs": [],
   "source": [
    "# padding\n",
    "def pad_with(vector, pad_width, iaxis, kwargs):\n",
    "    pad_value = kwargs.get('padder', 10.)\n",
    "    vector[:pad_width[0]] = pad_value\n",
    "    vector[-pad_width[1]:] = pad_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be6347a0-ec95-43ab-922c-b4791014bcae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cropping\n",
    "def crop_at(img,cropx,cropy,off7x=0,off7y=0):\n",
    "    y,x = img.shape\n",
    "    startx = x//2-(cropx//2)+off7x\n",
    "    starty = y//2-(cropy//2)+off7y\n",
    "    return img[starty:starty+cropy,startx:startx+cropx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c091213-6a7b-487f-ac42-92dd16dcabd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#clean too low flaot : here to delete imaginary part\n",
    "def clean_tiny_cplx(x, tol=1e-12):\n",
    "  tol = tol/ np.finfo(float).eps\n",
    "  x = np.real_if_close(x,tol=tol)\n",
    "  return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd0f05f5-8ca3-4241-b502-ecd1e38982d3",
   "metadata": {},
   "source": [
    "# Make a bullseye 32x32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4a2ce1e-d509-4164-b127-f354ce3c56e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nin = 32\n",
    "f = 2./Nin\n",
    "imgOrig = np.zeros((Nin,Nin),dtype=float)\n",
    "def trig(x):\n",
    "    mask1 = 3*x<1.\n",
    "    mask2 = 3*x>=2.\n",
    "    mask3 = 3*x<3.\n",
    "    if mask1|(mask2&mask3):\n",
    "        return 1.\n",
    "    else:\n",
    "        return 0.\n",
    "for iy in range(-Nin//2,Nin//2):\n",
    "    for ix in range(-Nin//2,Nin//2):\n",
    "        fx = f*ix\n",
    "        fy = f*iy\n",
    "        rad = np.sqrt(fx**2+fy**2)\n",
    "        imgOrig[ix+15,iy+15]=trig(rad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d31e5a20-15c3-4208-921a-f2616e3cf804",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotImg(img,fig, title=None, vmin=None, vmax=None):\n",
    "    ax = axes()\n",
    "    im = imshow(img, cmap='gray_r',vmin=vmin, vmax = vmax)\n",
    "    plt.title(title)\n",
    "    # Create new axes according to image position\n",
    "    cax = fig.add_axes([ax.get_position().x1+0.01,\n",
    "                    ax.get_position().y0,\n",
    "                    0.05,\n",
    "                    ax.get_position().height])\n",
    " \n",
    "    # Plot vertical colorbar\n",
    "    colorbar(im, cax=cax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d41d747a-d270-4d92-bef2-0cec60eefc0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(4,4))\n",
    "plotImg(imgOrig,fig, title=\"Original image 32x32\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a85d49a-39c4-4d19-85cf-d72e04ab5308",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(4,4))\n",
    "plotImg(imgOrig2,fig, title=\"Original image 32x32\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c90e0c50-f8ac-4e56-a5c6-53c99d9e0120",
   "metadata": {},
   "source": [
    "#  Pad to get a 128x128 image before x->u DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9333659-2541-43a5-a5ec-3731362249be",
   "metadata": {},
   "outputs": [],
   "source": [
    "orig_pad_128 = np.pad(imgOrig, (128-32)//2, pad_with, padder=0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3be48875-5cad-4305-8e0b-975cc506bd89",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(orig_pad_128,fig,title=\"Original image padded to 128x128\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "023f2dd1-690a-4cd2-b074-126f7bec6dbc",
   "metadata": {},
   "source": [
    "# Paper Eq 9 in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7e4e121-06d9-4252-a273-a8959f82baa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# specific for the periodic cardinal sinus (the best)\n",
    "def fhat2DSincWrapped(img,u1,u2):\n",
    "  # apply DFT to img\n",
    "  N1,N2 = img.shape\n",
    "  atildekl = npdft2Dshift(img) # formule (6) of fourier coeff in 2D\n",
    "  mtx1 = np.array([KuSincWrapped(u-k1/N1,N1) for k1 in range(-N1//2,N1//2) for u in u1]).reshape(N1,-1)\n",
    "  mtx2 = np.array([KuSincWrapped(v-k2/N2,N2) for k2 in range(-N2//2,N2//2) for v in u2]).reshape(N2,-1)\n",
    "  #print(mtx1.shape,fkl. shape,mtx2.shape)\n",
    "  return mtx1.T @ atildekl @ mtx2    # Formule (9) in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4808f94c-b0f4-4d16-8ae2-32247f28ca1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generic case\n",
    "def fhatKernel(img,u1,u2,kerWrapped):\n",
    "  # apply DFT to img\n",
    "  N1,N2 = img.shape\n",
    "  atildekl = npdft2Dshift(img) # formule (6) of fourier coeff in 2D\n",
    "  mtx1 = np.array([kerWrapped(u-k1/N1,N1) for k1 in range(-N1//2,N1//2) for u in u1]).reshape(N1,-1)\n",
    "  mtx2 = np.array([kerWrapped(v-k2/N2,N2) for k2 in range(-N2//2,N2//2) for v in u2]).reshape(N2,-1)\n",
    "  #print(mtx1.shape,fkl. shape,mtx2.shape)\n",
    "  return mtx1.T @ atildekl @ mtx2    # Formule (9) in 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2678e253-ec84-497d-878b-836cbabc22f7",
   "metadata": {},
   "source": [
    "# Define a new u-grid for Fourier resampling "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd03ac3b-db35-4ebc-a98c-ff286d5882dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "N1pad,N2pad = orig_pad_128.shape\n",
    "assert N1pad==N2pad, \"we suppose after that height = width\"\n",
    "Npad = N1pad\n",
    "Npad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f51c6a04-414c-4c8a-8035-75b70daf76c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "M=Npad*4 #width/height of the final image (Figure 2 if bulleteyes is 32x32)\n",
    "# If M=Npad/n n >= 1 the ghost images pile-up on the central image \n",
    "dx=1\n",
    "du = 1/(M*dx)\n",
    "umax = (M/2)*du\n",
    "umax,du"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45f70373-a7f0-4188-b4f0-f7f6ec980147",
   "metadata": {},
   "outputs": [],
   "source": [
    "u1=np.arange(-umax,umax,du)\n",
    "u2=np.arange(-umax,umax,du)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cab3ee65-b24e-49d0-9f4a-4952ab8890cd",
   "metadata": {},
   "source": [
    "# TF[Kx] kernel matrix in 2D Fourier grid\n",
    "$K_x(x)$ is the  Lanczos 3rd order renormalized (see interp1DKernels.py). \n",
    "$\\hat{K}_x(u)$ notation for the Fourier transform "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bf2bb30-8f46-4e25-9e5e-850953f92ace",
   "metadata": {},
   "outputs": [],
   "source": [
    "#generic\n",
    "def Kxhat2D(u1,u2,ker):\n",
    "  # warning interpolation regular integer grid\n",
    "  return np.array([ker(u)*ker(v) for u in u1 for v in u2]).reshape(u1.shape[0],u2.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b369c42-5e9e-4333-acbe-ac5f5ae5da8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a little as this is not optimized\n",
    "Kxhat_mtx = Kxhat2D(u1,u2,hatlanczos3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01d74ad5-8616-4c82-9774-1a033aa78c49",
   "metadata": {},
   "source": [
    "# What is the result if  Ku is the \"ideal\" Sinc-periodic kernel?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c87d69a5-275f-420d-a3f0-284443f3313a",
   "metadata": {},
   "source": [
    "## Eqs 7-9 : Fourier Transform of $f_d(x)$\n",
    "$$\n",
    "\\hat{f}_d(u) = \\sum_{k=-N/2}^{N/2-1} \\hat{a}_k  K_u(u-k/N)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b3973ff-948a-42b6-b92f-108f8d32fac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_mtx_ideal = fhat2DSincWrapped(orig_pad_128,u1,u2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e930604a-36fe-4ce2-a1d3-f28bbbd11ba1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(np.abs(fhat_mtx_ideal),fig,title=\"$\\hat{f}_d(u)$ (ideal $K_u$)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b62288a3-43a2-41cf-8038-94a2eb867b9c",
   "metadata": {},
   "source": [
    "## Eq. 11: convolution in Fourier space $\\hat{F}(u) = \\hat{K}_x(u) \\times \\hat{f}_d(u)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff1c5ed6-cd80-43e3-b6c4-a4558a58b352",
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_mtx_ideal = Kxhat_mtx * fhat_mtx_ideal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "790f43cc-5115-4b87-9794-bd405b37bbce",
   "metadata": {},
   "source": [
    "## $F(x)$ as inverse transform of $\\hat{F}(u)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "242755b8-8179-4304-a18c-3ffe497eabf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_real_ideal = npinvdft2Dshift(Fhat_mtx_ideal)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24be536f-12f8-41ee-bf16-badb621353ed",
   "metadata": {},
   "source": [
    "## cleaning for small numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df3a42f6-8c3f-4e45-85b8-fc0e739dfdf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_real_ideal = clean_tiny_cplx(F_real_ideal)\n",
    "F_real_ideal[np.abs(F_real_ideal)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b6e9ea3-d171-4ead-997b-7fb8501d247c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(F_real_ideal,fig,title=r\"$F(x)$ with $K_u=$Sinc-periodic\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7e6c974-dbc8-4ae0-8208-9c48a5a59b1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(F_real_ideal,fig,title=r\"$F(x)$ with $K_u=$Sinc-periodic\",\n",
    "        vmin=-0.003,vmax=0.003)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e093e62d-6480-47ab-8bd8-22c1c1db5d11",
   "metadata": {},
   "source": [
    "# Replace ideal $K_u$ kernel by lanczos3 (renormalized) periodic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71b09056-4146-4fa2-90c3-5238661e3605",
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_lanczos3_mtx = fhatKernel(orig_pad_128,u1,u2,lanczosRenorm3Wrapped)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "729b0183-85ce-4c91-86e6-92b03c611a4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_lanczos3_mtx = Kxhat_mtx * fhat_lanczos3_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "173ff798-ead7-424b-867a-898072b9363b",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_lanczos3_real = npinvdft2Dshift(Fhat_lanczos3_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64816704-77f3-40c0-8728-ba68436c5af0",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_lanczos3_real = clean_tiny_cplx(F_lanczos3_real)\n",
    "F_lanczos3_real[np.abs(F_lanczos3_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec89dbf8-a54c-4042-a577-3b44f89b9212",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(F_lanczos3_real,fig,title=r\"$F(x)$ with $K_u=$Lanczos3\",\n",
    "        vmin=-0.003,vmax=0.003)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "439d0d2c-195d-4d6e-a790-5969ec2f116e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_lanczos3_real,64,64),fig,\n",
    "    title=\"Central image ($K_u=$Lanczos3)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30356da3-d3c8-42d2-abc9-ab1d5f6dea67",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_lanczos3_real,64,64,128,0),fig,vmin=-0.003,vmax=0.003,\n",
    "    title=\"ghost image (right) ($K_u=$Lanczos3)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6afbc72d-4320-4557-8c20-d43d322d85a9",
   "metadata": {},
   "source": [
    "# Replace ideal $K_u$ kernel by Cubic periodic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93fecab8-896d-44f4-9d22-72dd607a2e77",
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_cubic_mtx = fhatKernel(orig_pad_128,u1,u2,cubicWrapped) # Ku = cubic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e14f0913-e6bd-4b77-bbd4-f4c11fd0f8e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_cubic_mtx = Kxhat_mtx * fhat_cubic_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98770bcd-4cd7-460f-b8aa-f688428a0d36",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_cubic_real = npinvdft2Dshift(Fhat_cubic_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7d37093-da95-4946-8b23-25b459fab928",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_cubic_real = clean_tiny_cplx(F_cubic_real)\n",
    "F_cubic_real[np.abs(F_cubic_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4baccbc-49fa-410e-a6e8-d31596039dfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(F_cubic_real,fig,title=r\"$F(x)$ with $K_u=$Cubic\",\n",
    "        vmin=-0.003,vmax=0.003)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae1e9b11-8bbf-4b33-ac8b-886e7f7279c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_cubic_real,64,64),fig,\n",
    "    title=\"Central image ($K_u=$Cubic)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfdf929a-4fc1-492c-891c-151a87a6ff93",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_cubic_real,64,64,128,0),fig,vmin=-0.003,vmax=0.003,\n",
    "    title=\"ghost image (right) ($K_u=$Cubic)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "128d07de-bbd1-42b6-85c0-de0e0ba80541",
   "metadata": {},
   "source": [
    "# Replace ideal $K_u$ kernel by Quintic periodic (Bernstein & Gruen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75684b84-e4d3-4bde-8836-de1bc92b7b0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_quinticBG_mtx = fhatKernel(orig_pad_128,u1,u2,quinticBGWrapped) # Ku = Quintic BG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f0b2c5c-a7b4-4d87-8d60-ab330e25f6c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_quinticBG_mtx = Kxhat_mtx * fhat_quinticBG_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ddd989f-65b4-435b-a742-b097aba22c66",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticBG_real = npinvdft2Dshift(Fhat_quinticBG_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fac881aa-4f5d-4904-b2a1-3c76ee6c866a",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticBG_real = clean_tiny_cplx(F_quinticBG_real)\n",
    "F_quinticBG_real[np.abs(F_quinticBG_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79cf91cd-8523-432c-9a86-d5a53e02abce",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(F_quinticBG_real,fig,title=r\"$F(x)$ with $K_u=$Quintic (BG)\",\n",
    "        vmin=-0.001,vmax=0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aedd23c-05ee-4c7b-9ea7-aaebc3f4cb44",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticBG_real,64,64),fig,\n",
    "    title=\"Central image ($K_u=$Quintic (BG))\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afab2f37-6b21-40b1-a492-e70901f6dcdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticBG_real,64,64,128,0),fig,vmin=-0.001,vmax=0.001,\n",
    "    title=\"ghost image (right) ($K_u=$Quintic (BG))\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76d30f61-f04f-4b7e-a727-cc9b68901181",
   "metadata": {},
   "source": [
    "# Replace ideal $K_u$ kernel by Quintic periodic (JE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30fa8705-577d-44de-a55a-29d3c7ca93cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_quinticJE_mtx = fhatKernel(orig_pad_128,u1,u2,quinticJEWrapped) # Ku = Quintic JE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf70c69b-d0c8-4ee5-9084-270a9202ca23",
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_quinticJE_mtx = Kxhat_mtx * fhat_quinticJE_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e3e38d9-ae54-42e4-9c58-31eaa60deca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticJE_real = npinvdft2Dshift(Fhat_quinticJE_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac4d3478-a798-4fa7-9168-e22c49591e94",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticJE_real = clean_tiny_cplx(F_quinticJE_real)\n",
    "F_quinticJE_real[np.abs(F_quinticJE_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3292e7a0-86c7-4d21-ba72-3bd78192e153",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(5,5))\n",
    "plotImg(F_quinticJE_real,fig,title=r\"$F(x)$ with $K_u=$Quintic (JE)\",\n",
    "        vmin=-0.001,vmax=0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d45bb655-9cdb-40d4-b8ac-5c321de9935a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticJE_real,64,64),fig,\n",
    "    title=\"Central image ($K_u=$Quintic (JE))\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c455f425-06c3-4c9a-91f8-a3d832e375a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticJE_real,64,64,128,0),fig,vmin=-0.001,vmax=0.001,\n",
    "    title=\"ghost image (right) ($K_u=$Quintic (JE))\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c4d4cc6-c6ef-4889-8243-f07de1b0936c",
   "metadata": {},
   "source": [
    "# F\\<Kernel\\> - F\\<Sinc Periodic\\> in real space (all and closeup on the central and 1st ghost right most)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68ffeb69-132d-4476-91c8-7dc9a48b45c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# common LUT\n",
    "vmin = -0.003 #-0.003\n",
    "vmax = 0.003  #0.003"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9a12387-a033-4b45-a686-84e9bd8ee5a8",
   "metadata": {},
   "source": [
    "## Lanczos-3 norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24d0baed-e454-4f62-9f0f-8f74d1bf4cc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(5,5))\n",
    "plotImg(F_lanczos3_real-F_real_ideal, fig,vmin=vmin, vmax=vmax, title=\"Lanczos-3 norm - Ideal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab872cfd-88e7-429a-aa42-6c23eb3f1d29",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_lanczos3_real-F_real_ideal,64,64),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Lanczos-3 norm - Ideal central\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33b413bb-b674-405e-8752-ba664dd4cf34",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_lanczos3_real-F_real_ideal,64,64,128,0),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Lanczos-3 norm - Ideal 1st ghost\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b0616df-d60a-4c7a-8803-a1c03469f9c2",
   "metadata": {},
   "source": [
    "## Cubic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "427894cf-689e-46fe-bc66-89f406b90a81",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(5,5))\n",
    "plotImg(F_cubic_real-F_real_ideal, fig,vmin=vmin, vmax=vmax, title=\"Cubic - Ideal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10a038f0-e19e-4d99-b9d1-722830b7f3ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_cubic_real-F_real_ideal,64,64),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Cubic - Ideal central\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac46b66b-5e9c-47ec-aae7-20d8a0b06a1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_cubic_real-F_real_ideal,64,64,128,0),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Cubic - Ideal 1st ghost\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a324f4-0671-45c0-ac7f-b608e9b2406f",
   "metadata": {},
   "source": [
    "## Quintic  (Bernstein & Gruen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a24e0341-c644-4ae4-8ce3-ebe4e2c62743",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(5,5))\n",
    "plotImg(F_quinticBG_real-F_real_ideal, fig,vmin=vmin, vmax=vmax, title=\"Quintic (BG) - Ideal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68d1e535-9440-4070-8432-7982092515e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticBG_real-F_real_ideal,64,64),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Quintic (BG) - Ideal central\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4e163cd-1532-42b7-95d3-749fa49b874e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticBG_real-F_real_ideal,64,64),fig,vmin=-0.001, vmax=0.001,\\\n",
    "    title=\"Quintic (BG) - Ideal central\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ba79e00-50ab-4650-9414-527a9db44e51",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticBG_real-F_real_ideal,64,64,128,0),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Quintic (BG) - Ideal 1st ghost\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dca16c4f-cce2-4253-8183-32f7a7f082ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticBG_real-F_real_ideal,64,64,128,0),fig,vmin=-0.001, vmax=0.001,\\\n",
    "    title=\"Quintic (BG) - Ideal 1st ghost\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d44a54c-f06f-4931-9cb8-cdc7fd051eeb",
   "metadata": {},
   "source": [
    "## Quintic  (JE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d34222-effa-4b6d-8f00-d17559e3dd89",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(5,5))\n",
    "plotImg(F_quinticJE_real-F_real_ideal, fig,vmin=vmin, vmax=vmax, title=\"Quintic (JE) - Ideal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3e6e98d-75b1-47a3-aee5-fc0f8297bde5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticJE_real-F_real_ideal,64,64),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Quintic (JE) - Ideal central\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d750607-abaf-4a47-beea-df1ae2ec358a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticJE_real-F_real_ideal,64,64),fig,vmin=-0.001, vmax=0.001,\\\n",
    "    title=\"Quintic (JE) - Ideal central\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f4ec16c-ca6e-463c-805d-b42aa48f3cf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticJE_real-F_real_ideal,64,64,128,0),fig,vmin=vmin, vmax=vmax,\\\n",
    "    title=\"Quintic (JE) - Ideal 1st ghost\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e62f6c4b-1f26-4160-baf0-4961c9bf8e86",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=figure(figsize=(4,4))\n",
    "plotImg(crop_at(F_quinticJE_real-F_real_ideal,64,64,128,0),fig,vmin=-0.001, vmax=0.001,\\\n",
    "    title=\"Quintic (JE) - Ideal 1st ghost\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "718d86a2-bef1-4001-b3d6-bedff7ed004b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}