image_interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/jecampagne/interpolation/blob/main/image_interpolation.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Xhujs6zr8BXM",
    "outputId": "86b5ccd7-d029-48ee-962d-a17336578fcd"
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "EQW8eI2CwtjG"
   },
   "outputs": [],
   "source": [
    "import matplotlib.gridspec as gridspec\n",
    "params = {'legend.fontsize':14,\n",
    "          'figure.figsize': (8, 8),\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,
   "metadata": {
    "id": "jq7i9HQOFYOE"
   },
   "outputs": [],
   "source": [
    "import matplotlib.ticker as mtick\n",
    "def decibel_formatter(unit='dB'):\n",
    "    def format_value(value, pos=None):\n",
    "        dB = 10*np.log10(value)\n",
    "        return f'{dB:.0f} {unit}'\n",
    "    return format_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "h5I8AefYxuhT"
   },
   "source": [
    "# Understanding the  Gary M. Bernstein & Daniel Gruen 2014 article\n",
    " see https://arxiv.org/pdf/1401.2636.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JJctpj_GFOpd"
   },
   "outputs": [],
   "source": [
    "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,
   "metadata": {
    "id": "O8wJtQXwFRKY"
   },
   "outputs": [],
   "source": [
    "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,
   "metadata": {
    "id": "KiKit9bvHQfR"
   },
   "outputs": [],
   "source": [
    "def fullprint(*args, **kwargs):\n",
    "  from pprint import pprint\n",
    "  import numpy\n",
    "  opt = numpy.get_printoptions()\n",
    "  numpy.set_printoptions(threshold=numpy.inf)\n",
    "  pprint(*args, **kwargs)\n",
    "  numpy.set_printoptions(**opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cQ1o-rSU5yBT"
   },
   "outputs": [],
   "source": [
    "def zero_tiny(x, threshold=1e-10):\n",
    "  x[np.abs(x) < threshold] = 0\n",
    "  return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "98_82bfDos1a"
   },
   "source": [
    "# Different interpolants in real & Fourier space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cMhy35jaqEgv"
   },
   "outputs": [],
   "source": [
    "def Unity(x,ker,Jmax=5):\n",
    "  js = np.arange(-Jmax,Jmax+1)\n",
    "  X,Y = np.meshgrid(x,js)\n",
    "  kxpj = ker(X+Y)\n",
    "  return np.sum(kxpj,axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Q4hOl73DW5bS"
   },
   "outputs": [],
   "source": [
    "# Whittaker-Shannon sinc. interpolant\n",
    "def sincInterp(x):\n",
    "  return np.sinc(x)\n",
    "def hatSincInterp(u):\n",
    "  u = np.abs(u)\n",
    "  return  np.piecewise(u,[u<0.5,u==0.5,u>0.5],[lambda x: 1., lambda x:0.5, lambda x:0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8MKrpG0qHQko"
   },
   "outputs": [],
   "source": [
    "def Bspline3(x):\n",
    "  # nb. degre poly is 3\n",
    "  x = np.abs(x)\n",
    "  def p1(x):\n",
    "    return 2/3+1/2* x * x * (x-2)\n",
    "  def p2(x):\n",
    "    return -1/6 *(x-2)*(x-2)*(x-2)\n",
    "  return np.piecewise(x, [x<1., (x>=1)&(x<2.), x>=2],[lambda x:p1(x),lambda x:p2(x),lambda x: 0.])\n",
    "\n",
    "def hatBspline3(u):\n",
    "  s = np.sinc(u)\n",
    "  s2 = s*s\n",
    "  return s2 * s2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Lwtl4uLQiPr8"
   },
   "outputs": [],
   "source": [
    "def hatBsplineCard3(u):\n",
    "  sinpu = np.sin(np.pi*u)\n",
    "  return hatBspline3(u)/(1-2/3*sinpu*sinpu)\n",
    "def s3(n):\n",
    "    n = n.astype(int)\n",
    "    r3 = np.sqrt(3)\n",
    "    return (-2*(n%2)+1)*r3/((2+r3)**(np.abs(n)))\n",
    "\n",
    "def BsplineCard3(xin):\n",
    "  xin = np.atleast_1d(xin)\n",
    "  xshape=xin.shape\n",
    "  x = xin.flatten()\n",
    "  nmin = np.ceil(x-2).astype(int)\n",
    "  nmax = np.floor(x+2).astype(int)\n",
    "  res = []\n",
    "  for i,x0 in enumerate(x):\n",
    "    res.append(np.sum(np.array([s3(n)*Bspline3(x0-n) for n in np.arange(nmin[i],nmax[i]+1,1)])))\n",
    "  return np.array(res).reshape(xshape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hatBspline5(u):\n",
    "  s = np.sinc(u)\n",
    "  s2 = s*s\n",
    "  return s2 * s2 * s2\n",
    "\n",
    "def hatBsplineCard5(u):\n",
    "  b1 = 13/60\n",
    "  b2 = 1/120\n",
    "  b0 = 1 - b1 - b2\n",
    "  piu = np.pi * u\n",
    "  return hatBspline5(u)/(b0 + 2*b1*np.cos(2*piu) +2*b2*np.cos(4*piu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 824
    },
    "id": "nJKXKJu6JUsi",
    "outputId": "80b35865-6bf7-4b40-d37e-a04f8cef6d37"
   },
   "outputs": [],
   "source": [
    "fig,axs = subplots(3,1,figsize=(8,8),tight_layout=True)\n",
    "\n",
    "kerX = [sincInterp,Bspline3,BsplineCard3]\n",
    "lab = [\"Sinc\",r\"$\\beta_3$\",r\"$\\gamma_3$\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "kerU = [hatSincInterp,hatBspline3,hatBsplineCard3]\n",
    "\n",
    "\n",
    "xv = np.linspace(-4,4,10_000)\n",
    "uvsmall = np.linspace(-2,2,10_000)\n",
    "uvlarge = np.linspace(0,5,10_000)\n",
    "#uvtiny = np.linspace(-0.1,0.1,10_000)\n",
    "for ik,(kX,kU) in enumerate(zip(kerX,kerU)):\n",
    "  axs[0].plot(xv,kX(xv),label=lab[ik],ls=lsStyle[ik])\n",
    "  axs[1].plot(uvsmall,kU(uvsmall),label=lab[ik],ls=lsStyle[ik])\n",
    "  axs[2].plot(uvlarge,np.abs(kU(uvlarge)),label=lab[ik],ls=lsStyle[ik])\n",
    "#  axs[3].plot(uvtiny,kU(1+uvtiny),label=lab[ik],ls=lsStyle[ik])\n",
    "\n",
    "for ax in axs:\n",
    "  ax.grid()\n",
    "axs[0].legend(loc='center right', bbox_to_anchor=(1.35, 0.5));\n",
    "#axs[1].legend()\n",
    "axs[1].set_xlabel(\"u\")\n",
    "\n",
    "axs[2].set_yscale(\"log\")\n",
    "axs[2].set_ylim([1e-7,2]);\n",
    "axs[2].set_xlabel(\"u\")\n",
    "axs[2].yaxis.set_major_formatter(mtick.FuncFormatter(decibel_formatter()))\n",
    "\n",
    "\n",
    "axs[0].set_xlabel(\"x\")\n",
    "axs[0].axhline(0., linestyle='--', color='black')\n",
    "axs[0].axvline(0., linestyle='--', color='black')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gjDYRKJTMn1b"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 400
    },
    "id": "Wn82wOQcPG5M",
    "outputId": "8ebcd520-abe6-4fd9-b28c-e04563091e63"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "xs=np.linspace(-3,3,1000,endpoint=True)\n",
    "beta4U = Unity(xs,Bspline3)\n",
    "plot(xs,beta4U,label=r\"$\\sum_j \\beta_3(x+j)$\")\n",
    "for j in np.arange(-5,6):\n",
    "  plot(xs,Bspline3(xs+j))#,label=fr\"$ \\beta_4(x+{j}) $\")\n",
    "grid()\n",
    "xlabel(\"x\")\n",
    "legend(loc=\"upper left\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 400
    },
    "id": "SqS9uRIv1dqu",
    "outputId": "f1ec4f5e-5e1b-42fc-c5ca-af4e138583d0"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "xs=np.linspace(-3,3,1000,endpoint=True)\n",
    "beta4U = Unity(xs,BsplineCard3)\n",
    "plot(xs,beta4U,label=r\"$\\sum_j \\gamma_3(x+j)$\")\n",
    "for j in np.arange(-5,6):\n",
    "  plot(xs,BsplineCard3(xs+j))#,label=fr\"$ \\gamma_4(x+{j}) $\")\n",
    "grid()\n",
    "xlabel(\"x\")\n",
    "#xlim([-3,3])\n",
    "legend(loc=\"upper left\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 400
    },
    "id": "p6SOB5QPMoKx",
    "outputId": "866eba65-4643-4225-dc82-b2f15cf3ae23"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "xs=np.linspace(-3,3,1000,endpoint=True)\n",
    "hatbeta4U = Unity(xs,hatBsplineCard3)\n",
    "plot(xs,hatbeta4U,label=r\"$\\sum_j \\hat{\\gamma}_3(x+j)$\")\n",
    "for j in np.arange(-5,6):\n",
    "  plot(xs,hatBspline3(xs+j))#,label=fr\"$ \\gamma_4(x+{j}) $\")\n",
    "grid()\n",
    "xlabel(\"x\")\n",
    "legend(loc=\"upper left\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Lvf-NMM0BmMO"
   },
   "outputs": [],
   "source": [
    "from scipy.special import sici\n",
    "\n",
    "# Lanczos interpolant of degree m as well as the background conservation corrected version\n",
    "def box(x):\n",
    "  x = np.abs(x)\n",
    "  return np.piecewise(x,[x<0.5,x==0.5,x>0.5],[lambda x: 1., lambda x:0.5, lambda x:0])\n",
    "\n",
    "def siInt(x):\n",
    "  si, _ = sici(x)\n",
    "  return si\n",
    "\n",
    "def lanczos(x,m):\n",
    "  x = np.abs(x)\n",
    "  return np.sinc(x)*np.sinc(x/m)*box(x/(2*m))\n",
    "\n",
    "def hatlanczos(u,m):\n",
    "  # the formula A12 are not correct.\n",
    "  # Here Galsim code  https://github.com/GalSim-developers/GalSim/blob/releases/2.5/src/Interpolant.cpp\n",
    "  # and validated with Mathematica\n",
    "  u = np.abs(u)\n",
    "  vp = m*(2*u+1)\n",
    "  vm = m*(2*u-1)\n",
    "  arg1 = vp+1\n",
    "  arg2 = vp-1\n",
    "  arg3 = vm-1\n",
    "  arg4 = vm+1\n",
    "  res = arg1*siInt(np.pi*arg1) - arg2*siInt(np.pi*arg2) + arg3*siInt(np.pi*arg3) - arg4*siInt(np.pi*arg4)\n",
    "  return res/(2*np.pi)\n",
    "\n",
    "def lanczos_norm(x,m):\n",
    "  c1 = hatlanczos(1,m)\n",
    "  return lanczos(x,m)*(1-2*c1*(np.cos(2*np.pi*x)-1)) # formula (24)\n",
    "\n",
    "def hatlanczos_norm(u,m):\n",
    "  c1 = hatlanczos(1,m)\n",
    "  return (1+2*c1)*hatlanczos(u,m) - c1*(hatlanczos(u+1,m)+hatlanczos(u-1,m))  # formula (25)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1V-bkhSa9gi-"
   },
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "lanczos_raw3 = partial(lanczos,m=3)\n",
    "hatlanczos_raw3 = partial(hatlanczos,m=3)\n",
    "lanczos_raw4 = partial(lanczos,m=4)\n",
    "hatlanczos_raw4 = partial(hatlanczos,m=4)\n",
    "lanczos_raw5 = partial(lanczos,m=5)\n",
    "hatlanczos_raw5 = partial(hatlanczos,m=5)\n",
    "\n",
    "\n",
    "lanczos3 = partial(lanczos_norm,m=3)\n",
    "hatlanczos3 = partial(hatlanczos_norm,m=3)\n",
    "\n",
    "lanczos4 = partial(lanczos_norm,m=4)\n",
    "hatlanczos4 = partial(hatlanczos_norm,m=4)\n",
    "\n",
    "lanczos5 = partial(lanczos_norm,m=5)\n",
    "hatlanczos5 = partial(hatlanczos_norm,m=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 423
    },
    "id": "wDpUdh2StgQ7",
    "outputId": "4dbdabae-1dfd-42bd-878b-a7d5f9dda0c4"
   },
   "outputs": [],
   "source": [
    "fig = figure(figsize=(6,4))\n",
    "ivals = np.arange(1,5)\n",
    "kerU = [hatlanczos_raw3,hatlanczos_raw5]\n",
    "lab = [\"Lanczos-3\",\"Lanczos-5\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "for  ik,ker in enumerate(kerU):\n",
    "  scatter(ivals,ker(ivals),label=lab[ik],ls=lsStyle[ik])\n",
    "legend()\n",
    "grid()\n",
    "yscale(\"log\")\n",
    "xlabel(r\"$j$\")\n",
    "ylabel(r\"$\\hat{L}_m(j)$\")\n",
    "ylim([3e-7,3e-3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "t8kt29ug39vU",
    "outputId": "143964d6-e96f-4568-b841-76b5824f1a80"
   },
   "outputs": [],
   "source": [
    "[hatlanczos(0,m) for m in [3,4,5,6,7,8]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kdP3ekbIFdcJ"
   },
   "outputs": [],
   "source": [
    "def LanczosSum(x,n,m):\n",
    "  res = []\n",
    "  for val in x:\n",
    "    s = np.sum(lanczos(val+n,m))\n",
    "    res.append(s)\n",
    "  return np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "exdVDC-gFmSE"
   },
   "outputs": [],
   "source": [
    "xs = np.linspace(-5,5,1000,endpoint=True)\n",
    "data3 = LanczosSum(xs,np.arange(-5000,5001),3)\n",
    "data5 = LanczosSum(xs,np.arange(-5000,5001),5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 682
    },
    "id": "lI--VVuLGHOS",
    "outputId": "75fa96c9-5457-415e-a143-1b5470e82c26"
   },
   "outputs": [],
   "source": [
    "plot(xs,data3)\n",
    "plot(xs,data5)\n",
    "grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aYAtXol9GsV0"
   },
   "outputs": [],
   "source": [
    "def NormLanczosSum(x,n,m):\n",
    "  res = []\n",
    "  for val in x:\n",
    "    s = np.sum(lanczos_norm(val+n,m))\n",
    "    res.append(s)\n",
    "  return np.array(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "nEgh8LpHI9_N",
    "outputId": "4d6c3450-b425-492d-bb58-4278e5a0146d"
   },
   "outputs": [],
   "source": [
    "1- NormLanczosSum(np.array([0,0.25,0.5,0.75]),np.arange(-5000,5001),3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "VOmBo54tJOtZ",
    "outputId": "f142f075-92fd-4eac-a782-3dded7bb5312"
   },
   "outputs": [],
   "source": [
    "1-NormLanczosSum(np.array([0,0.25,0.5,0.75]),np.arange(-5000,5001),5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kTibZxhPG0za"
   },
   "outputs": [],
   "source": [
    "xs = np.linspace(-5,5,1000,endpoint=True)\n",
    "data3n = NormLanczosSum(xs,np.arange(-5000,5001),3)\n",
    "data5n = NormLanczosSum(xs,np.arange(-5000,5001),5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 477
    },
    "id": "5ZQxZl5WHGQW",
    "outputId": "75254ce6-371e-4786-99a2-59b5276f79a7"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(7,5))\n",
    "plot(xs,data3,label=\"Lanczos 3\")\n",
    "plot(xs,data5,label=\"Lanczos 5\")\n",
    "plot(xs,data3n,label=\"Lanczos Norm 3\")\n",
    "plot(xs,data5n,label=\"Lanczos Norm 5\")\n",
    "grid();\n",
    "legend()\n",
    "xlabel(r\"$x$\")\n",
    "ylabel(r\"$\\sum_{j=-\\infty}^{\\infty} K(x-j)$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_2n2dSReihU8"
   },
   "outputs": [],
   "source": [
    "def hatLanczosSum(u,n,m):\n",
    "  res = []\n",
    "  for val in u:\n",
    "    s = np.sum(hatlanczos(val+n,m))\n",
    "    res.append(s)\n",
    "  return np.array(res)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7PZMP0B10ZmC"
   },
   "outputs": [],
   "source": [
    "us = np.linspace(-0.5,0.5,endpoint=True)\n",
    "data3 = hatLanczosSum(us,np.arange(-3000,3001),3)\n",
    "data5 = hatLanczosSum(us,np.arange(-3000,3001),5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 721
    },
    "id": "aW8fz2PR2PQb",
    "outputId": "93b58a7d-be2d-4d54-e98a-45e97cea3062"
   },
   "outputs": [],
   "source": [
    "plot(us,data3-1)\n",
    "plot(us,data5-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "gx1-HWuLFdAJ",
    "outputId": "f00f6cad-f6ff-4fc1-e4f2-e3b7ebb2a343"
   },
   "outputs": [],
   "source": [
    "ivals = np.arange(-10,11)\n",
    "ivals,hatlanczos(ivals,3), hatlanczos3(ivals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bAsBv9OiekqZ"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nJKnLvErbxRG"
   },
   "outputs": [],
   "source": [
    "def PlotKer(kerX,kerU,lab,lsStyle,ax20ylim,ax21ylim):\n",
    "  fig = plt.figure(tight_layout=True,figsize=(8,10))\n",
    "  gs = gridspec.GridSpec(3, 2)\n",
    "  ax0 = fig.add_subplot(gs[0, :])\n",
    "  ax1 = fig.add_subplot(gs[1, :])\n",
    "  ax20 = fig.add_subplot(gs[2, 0])\n",
    "  ax21 = fig.add_subplot(gs[2, 1])\n",
    "  axs = [ax0,ax1,ax20,ax21]\n",
    "\n",
    "  xv = np.linspace(-4,4,10_000)\n",
    "  uvsmall = np.linspace(0,3,10_000)\n",
    "  uvlarge = np.linspace(0,5,10_000)\n",
    "  uvtiny1 = np.linspace(0.8,1.2,10_000)\n",
    "  uvtiny0 = np.linspace(0,0.5,10_000)\n",
    "\n",
    "  for ik,(kX,kU) in enumerate(zip(kerX,kerU)):\n",
    "    ax0.plot(xv,kX(xv),label=lab[ik],ls=lsStyle[ik])\n",
    "    ax1.plot(uvlarge,np.abs(kU(uvlarge)),label=lab[ik],ls=lsStyle[ik])\n",
    "    ax20.plot(uvtiny0,kU(uvtiny0),label=lab[ik],ls=lsStyle[ik])\n",
    "    ax21.plot(uvtiny1,kU(1+uvtiny1),label=lab[ik],ls=lsStyle[ik])\n",
    "\n",
    "  for ax in axs:\n",
    "    ax.grid()\n",
    "  ax0.legend()#loc='center right', bbox_to_anchor=(1.35, 0.5));\n",
    "  ax0.set_xlabel(\"x\")\n",
    "  ax0.axhline(0., linestyle='--', color='black')\n",
    "  ax0.axvline(0., linestyle='--', color='black')\n",
    "\n",
    "  ax1.set_yscale(\"log\")\n",
    "  ax1.set_xlabel(\"u\")\n",
    "  ax1.yaxis.set_major_formatter(mtick.FuncFormatter(decibel_formatter()))\n",
    "  ax1.set_ylim([1e-7,2]);\n",
    "  ax20.set_xlabel(\"u\")\n",
    "  ax20.set_ylim(ax20ylim);\n",
    "  ax21.set_xlabel(\"u\");\n",
    "  ax21.set_ylim(ax21ylim);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "kdOrup5hdHSl",
    "outputId": "b7ab2729-c472-411e-a83f-4fffb1259d68"
   },
   "outputs": [],
   "source": [
    "kerX = [sincInterp,lanczos_raw3,lanczos_raw4,lanczos_raw5]\n",
    "lab = [\"Sinc\",\"Lanczos-3\",\"Lanczos-4\",\"Lanczos-5\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "kerU = [hatSincInterp,hatlanczos_raw3,hatlanczos_raw4,hatlanczos_raw5]\n",
    "\n",
    "PlotKer(kerX,kerU,lab,lsStyle,ax20ylim=[0.98,1.02],ax21ylim=[-2e-4,2e-4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "K2l0_oeQ4mt0"
   },
   "outputs": [],
   "source": [
    "# linear interpolant\n",
    "def h1(x):\n",
    "  x = np.abs(x)\n",
    "  return np.piecewise(x,[x<1.,x>=1.],[lambda x: 1-x, lambda x:0])\n",
    "def hat1(u):\n",
    "  u = np.abs(u)\n",
    "  s = np.sinc(u)\n",
    "  return s*s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9k9qcRUsfAjj"
   },
   "outputs": [],
   "source": [
    "# cubic interpolant\n",
    "def h3(x):\n",
    "  # warning this is interpolator with perfect match at interger values\n",
    "  x = np.abs(x)\n",
    "  def p1(x): # x<1\n",
    "    return 1+x*x*(1.5*x-2.5)\n",
    "  def p2(x):\n",
    "    return -0.5*(x-1.0)*(x-2.0)*(x-2.0)\n",
    "  return np.piecewise(x, [x<1., (x>=1)&(x<2.), x>=2],[lambda x:p1(x),lambda x:p2(x),lambda x: 0.])\n",
    "\n",
    "def hat3(u):\n",
    "  # warning this is interpolator with perfect match at interger values\n",
    "  u = np.abs(u)\n",
    "  s = np.sinc(u)\n",
    "  c = np.cos(np.pi * u)\n",
    "  return s * s * s * (3.0 * s - 2.0 * c)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_maZUAk4tWLk"
   },
   "outputs": [],
   "source": [
    "# quitic interpolants\n",
    "# Gary original paper version\n",
    "def h5Gary(x):\n",
    "  x = np.abs(x)\n",
    "  def p1(x): # x<1\n",
    "    return 1.0 + x*x*x*(-95./12. + x*(23./2. + x*(-55./12.)))\n",
    "  def p2(x): # 1<=x<2\n",
    "    return (x-1.)*(x-2.)*(-23./4. + x*(29./2. + x*(-83./8. + x*(55./24.))))\n",
    "  def p3(x): # 2<=x<3\n",
    "    return (x-2.)*(x-3.)*(x-3.)*(-9./4. + x*(25./12. + x*(-11./24.)))\n",
    "  return np.piecewise(x, [x<1., (x>=1)&(x<2.), (x>=2)&(x<3.), x>=3.],[lambda x:p1(x),lambda x:p2(x),lambda x:p3(x),lambda x: 0.])\n",
    "\n",
    "def hat5Gary(u): #x-checked by Mathematica\n",
    "  u = np.abs(u)\n",
    "  s = np.sinc(u)\n",
    "  s2 = s*s\n",
    "  s4 = s2*s2\n",
    "  s5 = s*s4\n",
    "  piu = np.pi * u\n",
    "  c = np.cos(piu)\n",
    "  piu2 = piu*piu\n",
    "  return s5 * (2.0 * (-27.0 + piu2) * c + s * (55.0 -19.0 * piu2))\n",
    "\n",
    "# Galsim ALT_QUINTIC by Mike Jarvis (BUT MJ indicates that ALT_QUINTIC is not used)\n",
    "def h5MJ(x):\n",
    "  x = np.abs(x)\n",
    "  def p1(x): # x<1\n",
    "    return 1. + x*x*x*(-15./2. + x*(32./3. + x*(-25./6.)))\n",
    "  def p2(x): # 1<=x<2\n",
    "    return (x-1.)*(x-2.)*(-23./4. + x*(169./12. + x*(-39./4. + x*(25./12.))))\n",
    "  def p3(x): # 2<=x<3\n",
    "    return (x-2.)*(x-3.)*(x-3.)*(x-3.)*(3./4. + x*(-5./12.))\n",
    "  return np.piecewise(x, [x<1., (x>=1)&(x<2.), (x>=2)&(x<3.), x>=3.],[lambda x:p1(x),lambda x:p2(x),lambda x:p3(x),lambda x: 0.])\n",
    "\n",
    "def hat5MJ(u):\n",
    "  u = np.abs(u)\n",
    "  s = np.sinc(u)\n",
    "  s2 = s*s\n",
    "  s4 = s2*s2\n",
    "  piu = np.pi * u\n",
    "  c = np.cos(piu)\n",
    "  piu2 = piu*piu\n",
    "  return -s4*(s2*(12.*piu2-50.) + 44.*s*c + 5.);\n",
    "\n",
    "# Version by Erik H. W. Meijering, Karel J. Zuiderveld, and Max A. Viergever\n",
    "# IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 8, NO. 2, FEBRUARY 1999 p 192\n",
    "\n",
    "def h5MZVgen(x,c):\n",
    "  x = np.abs(x)\n",
    "  def p1(x): # x<1\n",
    "    return 1. + x*x*(-5./2.+8*c + x*x*(45./16.-18*c + x*(-21./6.+10.*c)))\n",
    "  def p2(x): # 1<=x<2\n",
    "    return 5. + x*(-15. + x*(35./2. + x*(-10. + (45./16. - (5.*x)/16.)*x))) + \\\n",
    "      c * (-66. + x*(265. + x*(-392. + x*(270. + x*(-88. + 11.*x)))))\n",
    "  def p3(x): # 2<=x<3\n",
    "    return c*(-162. + x*(297. + x*(-216. + x*(78. + (-14. + x)*x))))\n",
    "  return np.piecewise(x, [x<1., (x>=1)&(x<2.), (x>=2)&(x<3.), x>=3.],[lambda x:p1(x),lambda x:p2(x),lambda x:p3(x),lambda x: 0.])\n",
    "\n",
    "def h5MZV(x):\n",
    "  x = np.abs(x)\n",
    "  def p1(x): # x<1\n",
    "    return 1. + x*x * (-(17./8.) + (63./32. - (27.* x)/32.)*x*x)\n",
    "  def p2(x): # 1<=x<2\n",
    "    return 1./64.*(-2. + x) * (-1. + x) * (61. + x* (9. + x* (-45. + 13.* x)))\n",
    "  def p3(x): # 2<=x<3\n",
    "    return 3./64.* (-3. + x)*(-3. + x)*(-3. + x)*(-3. + x)*(-2. + x)\n",
    "  return np.piecewise(x, [x<1., (x>=1)&(x<2.), (x>=2)&(x<3.), x>=3.],[lambda x:p1(x),lambda x:p2(x),lambda x:p3(x),lambda x: 0.])\n",
    "\n",
    "\n",
    "def hat5MZV(u):\n",
    "  u = np.abs(u)\n",
    "  def fapp(u):\n",
    "    u2 = u*u\n",
    "    u4 = u2*u2\n",
    "    P = 0.33333333333333333333 - 0.24271404762585718748 * u2 + 0.039750660445573631826 * u4\n",
    "    Q = 1.0000000000000000000  + 0.25881829723136429944 * u2 + 0.026805791044710590408 * u4\n",
    "    return P/Q\n",
    "\n",
    "  def fno_app(u):\n",
    "    piu = np.pi * u\n",
    "    piu2 = piu*piu\n",
    "    s = np.sinc(u)\n",
    "    c = np.cos(piu)\n",
    "    return (s-c)/piu2\n",
    "\n",
    "  def f(u):\n",
    "    return np.piecewise(u,[u<0.1,u>=0.1],[lambda u: fapp(u), lambda u: fno_app(u)])\n",
    "\n",
    "  piu = np.pi * u\n",
    "  piu2 = piu*piu\n",
    "  c = np.cos(piu)\n",
    "  s = np.sinc(u)\n",
    "  s2 = s*s\n",
    "  s3 = s*s2\n",
    "  return 3./8. * s3 *( 35. * f(u) + s2 * (6.*c -15.*s) )\n",
    "\n",
    "# A variation upon Gary quintic interpolant\n",
    "# that I have cooked such that H^(l)[u=1]=0 for l in{1,...,5}\n",
    "# nb. Gary in contrast is defined by h^(2)[x=0]=0\n",
    "def h5JE(x):\n",
    "  x = np.abs(x)\n",
    "  def p1(x): # x<1\n",
    "    x2 = x * x\n",
    "    pi2 = np.pi * np.pi\n",
    "    return (15. *(-12. + x2*(27. + x*(-13. + (-3. + x)*x))) \\\n",
    "            + pi2 * (12. - x2*(15. + x*(35. + x*(-63. + 25.*x)))))/(12.*(-15. +pi2))\n",
    "  def p2(x): # 1<=x<2\n",
    "    pi2 = np.pi * np.pi\n",
    "    return ((-2. + x) * (-1 + x)*(-15.*(24. + x*(-3. + (-6. + x)*x))\\\n",
    "                + pi2 * (-48. + x * (153. + x*(-114. + 25.*x)))))/(24.*(-15. + pi2))\n",
    "  def p3(x): # 2<=x<3\n",
    "    pi2 = np.pi * np.pi\n",
    "    return -(((-3. + x)*(-3.+x)*(-2. + x)*(-3.* (-7. + x)* x\\\n",
    "                + pi2*(-3. + x)*(-8. + 5.*x)))/(24.*(-15. + pi2)))\n",
    "  return np.piecewise(x, [x<1., (x>=1)&(x<2.), (x>=2)&(x<3.), x>=3.],\\\n",
    "             [lambda x:p1(x),lambda x:p2(x),lambda x:p3(x),lambda x: 0.])\n",
    "\n",
    "def hat5JE(u):\n",
    "  u = np.abs(u)\n",
    "  pi2 = np.pi * np.pi\n",
    "  piu = np.pi * u\n",
    "  piu2 = piu*piu\n",
    "  c = np.cos(piu)\n",
    "  ss = np.sin(piu)\n",
    "  s = np.sinc(u)\n",
    "  s2 = s*s\n",
    "  s4 = s2*s2\n",
    "  s5 = s*s4\n",
    "  return (s5* (np.pi*(24.*np.pi* (-1. + u*u)*c - (39. + 7.*pi2)* u*ss)\\\n",
    "               + 5.*(-3. + 5.*pi2)*s))/(-15. + pi2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dOzM5eYcbFtA"
   },
   "outputs": [],
   "source": [
    "# Magic Kernel used in Facebook & Instagram https://johncostella.com/magic/mks.pdf\n",
    "def hatmagic3(u):\n",
    "  u = np.abs(u)\n",
    "  s = np.sinc(u)\n",
    "  piu = np.pi * u\n",
    "  sin = np.sin(piu)\n",
    "  return s*s*s/(1-0.5*sin*sin)\n",
    "def mks2021(x):\n",
    "  x = np.abs(x)\n",
    "  def p1(x): # x<1/2\n",
    "    return 577./576. - 239./144. * x*x\n",
    "  def p2(x): #1/2<=x<3/2\n",
    "    return (239. - 379. * x + 140. * x*x )/144.\n",
    "  def p3(x): #3/2<=x<5/2\n",
    "    return (-130. + 113. * x - 24. * x*x )/144.\n",
    "  def p4(x): #5/2<=x<7/2\n",
    "    return (45. - 27.* x + 4. * x*x )/144.\n",
    "  def p5(x): #7/2<=x<9/2\n",
    "    return (-81. + 36. * x - 4. * x*x)/1152.\n",
    "  return np.piecewise(x,\n",
    "            [x<0.5, (x>=0.5)&(x<1.5), (x>=1.5)&(x<2.5), (x>=2.5)&(x<3.5), (x>=3.5)&(x<4.5), x>=4.5],\n",
    "             [lambda x:p1(x),lambda x:p2(x),lambda x:p3(x),lambda x:p4(x),lambda x:p5(x),lambda x: 0.])\n",
    "\n",
    "def hatmks2021(u):\n",
    "  u = np.abs(u)\n",
    "  piu = np.pi *u\n",
    "  s = np.sinc(u)\n",
    "  return s*s*s*(102. - 35.*np.cos(2*piu) + 6.* np.cos(4*piu) - np.cos(6*piu))/72."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "R33ImNIhInqb",
    "outputId": "4ad7b165-34be-4885-aafb-b89dc0d1c7b0"
   },
   "outputs": [],
   "source": [
    "kerX = [sincInterp,h1,h3,h5MZV]\n",
    "lab = [\"Sinc\",\"Linear\",\"Cubic (MZV)\",\"Quintic (MZV)\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-1]=\"--\"\n",
    "kerU = [hatSincInterp,hat1,hat3,hat5MZV]\n",
    "\n",
    "PlotKer(kerX,kerU,lab,lsStyle,ax20ylim=[0.999,1.001],ax21ylim=[-2e-4,2e-4])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 676
    },
    "id": "744O_q0A8yey",
    "outputId": "624afffa-a91f-47aa-b15e-de18eacc409e"
   },
   "outputs": [],
   "source": [
    "fig,axs = subplots(3,1,figsize=(8,8))\n",
    "#kerX = [sincInterp,lanczos3,lanczos5,h1,h3,h5Gary,h5MZV,mks2021]\n",
    "#lab = [\"Sinc\",\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Linear\",\"Cubic\",\"Quintic (GB)\",\"Quintic (MZV)\",\"Quintic (JE)\",\"mks2021\"]\n",
    "#lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-1]=\"--\"\n",
    "#kerU = [hatSincInterp,hatlanczos3,hatlanczos5,hat1,hat3,hat5Gary,hat5MZV,hatmks2021]\n",
    "\n",
    "kerX = [sincInterp,lanczos3,lanczos5,h3,h5Gary,h5MZV,h5JE, mks2021]\n",
    "lab = [\"Sinc\",\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic\",\"Quintic (GB)\",\"Quintic (MZV)\",\"Quintic (JE)\",\"mks2021\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "lsStyle[-1]=\"--\"\n",
    "kerU = [hatSincInterp,hatlanczos3,hatlanczos5,hat3,hat5Gary,hat5MZV,hat5JE, hatmks2021]\n",
    "\n",
    "\n",
    "xv = np.linspace(-4,4,10_000)\n",
    "uvsmall = np.linspace(0,3,10_000)\n",
    "uvlarge = np.linspace(0,5,10_000)\n",
    "#uvtiny = np.linspace(-0.1,0.1,10_000)\n",
    "for ik,(kX,kU) in enumerate(zip(kerX,kerU)):\n",
    "  axs[0].plot(xv,kX(xv),label=lab[ik],ls=lsStyle[ik])\n",
    "  axs[1].plot(uvsmall,kU(uvsmall),label=lab[ik],ls=lsStyle[ik])\n",
    "  axs[2].plot(uvlarge,np.abs(kU(uvlarge)),label=lab[ik],ls=lsStyle[ik])\n",
    "#  axs[3].plot(uvtiny,kU(1+uvtiny),label=lab[ik],ls=lsStyle[ik])\n",
    "\n",
    "for ax in axs:\n",
    "  ax.grid()\n",
    "axs[0].legend(loc='center right', bbox_to_anchor=(1.35, 0.5));\n",
    "#axs[1].legend()\n",
    "axs[2].set_yscale(\"log\")\n",
    "axs[2].set_ylim([1e-7,2]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 474
    },
    "id": "qMYYif0Q4Awu",
    "outputId": "67648ee8-2e6a-4cd6-ccd0-86ddc0bdb627"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "\n",
    "lab = [\"Sinc\",\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic\",\"Quintic (MZV)\",\"Quintic (GB)\",\"Quintic (JE)\"]#, \"BSplineC 3\", \"BSplineC 5\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-2]=\"-.\"\n",
    "lsStyle[-1]=\"-.\"\n",
    "kerU = [hatSincInterp,hatlanczos3,hatlanczos5,hat3,hat5MZV,hat5Gary,hat5JE]\n",
    "\n",
    "#kerX = [sincInterp,lanczos3,lanczos5,h3,h5Gary,mks2021]\n",
    "#lab = [\"Sinc\",\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic\",\"Quintic (GB)\",\"mks2021\"]\n",
    "#lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-1]=\"--\"\n",
    "#kerU = [hatSincInterp,hatlanczos3,hatlanczos5,hat3,hat5Gary,hatmks2021]\n",
    "\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "\n",
    "uv = np.linspace(-0.3,0.3,1000)\n",
    "for ik,kU in enumerate(kerU):\n",
    "  plot(uv,kU(uv+1),label=lab[ik],ls=lsStyle[ik],c=colors[ik])\n",
    "ylim([-0.005,0.005])\n",
    "grid();\n",
    "legend(loc='center right', bbox_to_anchor=(1.45, 0.5),fontsize=12);\n",
    "# Fill in area under the curve and the horizontal lines\n",
    "fill_between(x=uv, y1=0.001, y2=-0.001, color='gray',  interpolate=True, alpha=.1)\n",
    "axvline(-0.25, linestyle='--', color='black')\n",
    "axvline(0.25, linestyle='--', color='black');\n",
    "axvline(-0.125, linestyle='--', color='black')\n",
    "axvline(0.125, linestyle='--', color='black');\n",
    "axvline(-0.5/6, linestyle='--', color='black')\n",
    "axvline(0.5/6, linestyle='--', color='black');\n",
    "ylabel(r\"$\\hat{K}(u+1)$\")\n",
    "xlabel(r\"$u=x/N$\");\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/6, -0.004), xycoords='data',\n",
    "    xytext=(0.5/6, -0.004), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x6 padding', xy=(0, -0.004), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points')\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/4, -0.003), xycoords='data',\n",
    "    xytext=(0.5/4, -0.003), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x4 padding', xy=(0, -0.003), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points')\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/2, -0.002), xycoords='data',\n",
    "    xytext=(0.5/2, -0.002), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x2 padding', xy=(0, -0.002), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 455
    },
    "id": "aLtd69rBCgbP",
    "outputId": "46596a56-a28c-4212-abb4-db72e59dbbea"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "\n",
    "lab = [\"Cubic\",\"Quintic (MZV)\", \"Quintic (GB)\",\"Quintic (JE)\",]#,\"BSplineC 3\", \"BSplineC 5\"]#,\"mks2021\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-2]=\"-.\"\n",
    "lsStyle[-1]=\"-.\"\n",
    "kerU = [hat3,hat5MZV,hat5Gary,hat5JE,]#,hatBsplineCard3,hatBsplineCard5]#,hatmks2021]\n",
    "\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "\n",
    "\n",
    "uv = np.linspace(-0.2,0.2,1000)\n",
    "for ik,kU in enumerate(kerU):\n",
    "  plot(uv,kU(uv+1),label=lab[ik],ls=lsStyle[ik],c=colors[ik+3])\n",
    "ylim([-0.003,0.003])\n",
    "grid();\n",
    "legend(loc='center right', bbox_to_anchor=(1.45, 0.5), fontsize=12);\n",
    "# Fill in area under the curve and the horizontal lines\n",
    "fill_between(x=uv, y1=0.0002, y2=-0.0002, color='blue',  interpolate=True, alpha=.1)\n",
    "fill_between(x=uv, y1=0.001, y2=-0.001, color='gray',  interpolate=True, alpha=.1)\n",
    "axvline(-0.125, linestyle='--', color='black')\n",
    "axvline(0.125, linestyle='--', color='black');\n",
    "#axvline(-0.25, linestyle='--', color='black')\n",
    "#axvline(0.25, linestyle='--', color='black');\n",
    "axvline(-0.5/6, linestyle='--', color='black')\n",
    "axvline(0.5/6, linestyle='--', color='black');\n",
    "ylabel(r\"$\\hat{K}(u+1)$\")\n",
    "xlabel(r\"$u=x/N$\");\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/6, -0.002), xycoords='data',\n",
    "    xytext=(0.5/6, -0.002), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x6 padding', xy=(0, -0.002), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points')\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/4, -0.0015), xycoords='data',\n",
    "    xytext=(0.5/4, -0.0015), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x4 padding', xy=(0, -0.0015), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points');\n",
    "\n",
    "#plt.annotate(\n",
    "#    '', xy=(-0.5/2, -0.0010), xycoords='data',\n",
    "#    xytext=(0.5/2, -0.0010), textcoords='data',\n",
    "#    arrowprops={'arrowstyle': '<->'})\n",
    "#plt.annotate(\n",
    "#    'x2 padding', xy=(0, -0.0010), xycoords='data', ha='center',\n",
    "#    xytext=(0, 1), textcoords='offset points')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 455
    },
    "id": "T_YbDxPlnV7b",
    "outputId": "7846b48c-0433-4bab-9f30-90027937ebd8"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "\n",
    "lab = [\"Sinc\",\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic\",\"Quintic (MZV)\", \"Quintic (GB)\",\"Quintic (JE)\"]#,\"BSplineC 3\",\"BSplineC 5\"]#\"mks2021\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-2]=\"-.\"\n",
    "lsStyle[-1]=\"-.\"\n",
    "kerU = [hatSincInterp,hatlanczos3,hatlanczos5,hat3,hat5MZV,hat5Gary,hat5JE]#,hatBsplineCard3, hatBsplineCard5]#hatmks2021]\n",
    "\n",
    "#kerX = [sincInterp,lanczos3,lanczos5,h3,h5Gary,mks2021]\n",
    "#lab = [\"Sinc\",\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic\",\"Quintic (GB)\",\"mks2021\"]\n",
    "#lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-1]=\"--\"\n",
    "#kerU = [hatSincInterp,hatlanczos3,hatlanczos5,hat3,hat5Gary,hatmks2021]\n",
    "\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "\n",
    "uv = np.linspace(-0.3,0.3,1000)\n",
    "for ik,kU in enumerate(kerU):\n",
    "  plot(uv,kU(uv+2),label=lab[ik],ls=lsStyle[ik],c=colors[ik])\n",
    "ylim([-0.01,0.01])\n",
    "grid();\n",
    "legend(loc='center right', bbox_to_anchor=(1.45, 0.5), fontsize=12);\n",
    "# Fill in area under the curve and the horizontal lines\n",
    "fill_between(x=uv, y1=0.001, y2=-0.001, color='gray',  interpolate=True, alpha=.1)\n",
    "axvline(-0.25, linestyle='--', color='black')\n",
    "axvline(0.25, linestyle='--', color='black');\n",
    "axvline(-0.125, linestyle='--', color='black')\n",
    "axvline(0.125, linestyle='--', color='black');\n",
    "axvline(-0.5/6, linestyle='--', color='black')\n",
    "axvline(0.5/6, linestyle='--', color='black');\n",
    "ylabel(r\"$\\hat{K}(u+2)$\")\n",
    "xlabel(r\"$u=x/N$\");\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/6, -0.009), xycoords='data',\n",
    "    xytext=(0.5/6, -0.009), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x6 padding', xy=(0, -0.009), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points')\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/4, -0.008), xycoords='data',\n",
    "    xytext=(0.5/4, -0.008), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x4 padding', xy=(0, -0.008), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points')\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(-0.5/2, -0.007), xycoords='data',\n",
    "    xytext=(0.5/2, -0.007), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x2 padding', xy=(0, -0.007), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qG_6WQHVGfUW"
   },
   "outputs": [],
   "source": [
    "def E0(u,ker,n=5):\n",
    "  sum=0\n",
    "  for j in range(1,n+1):\n",
    "    sum += ker(j+u)+ker(j-u)\n",
    "  return sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 455
    },
    "id": "Ts4FzBYoHJ1Q",
    "outputId": "b3e17169-0573-4b1c-dba8-238a6547c288"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "\n",
    "lab = [\"Lanczos Norm 5\", \"Cubic\",\"Quintic (MZV)\",\"Quintic (GB)\",\"Quintic (JE)\"]#,\"BSplineC 3\",\"BSplineC 5\"]#,\"mks2021\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-2]=\"-.\"\n",
    "lsStyle[-1]=\"-.\"\n",
    "kerU = [hatlanczos5,hat3,hat5MZV,hat5Gary,hat5JE]#,hatBsplineCard3, hatBsplineCard5]#,hatmks2021]\n",
    "\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "\n",
    "\n",
    "uv = np.linspace(0,0.2,1000)\n",
    "for ik,kU in enumerate(kerU):\n",
    "  plot(uv,E0(uv,kU),label=lab[ik],ls=lsStyle[ik],c=colors[ik+2])\n",
    "ylim([-0.002,0.002])\n",
    "grid();\n",
    "legend(loc='center right', bbox_to_anchor=(1.45, 0.5), fontsize=12);\n",
    "# Fill in area under the curve and the horizontal lines\n",
    "fill_between(x=uv, y1=0.0002, y2=-0.0002, color='blue',  interpolate=True, alpha=.1)\n",
    "fill_between(x=uv, y1=0.001, y2=-0.001, color='gray',  interpolate=True, alpha=.1)\n",
    "#axvline(-0.125, linestyle='--', color='black')\n",
    "axvline(0.125, linestyle='--', color='black');\n",
    "#axvline(-0.25, linestyle='--', color='black')\n",
    "#axvline(0.25, linestyle='--', color='black');\n",
    "#axvline(-0.5/6, linestyle='--', color='black')\n",
    "axvline(0.5/6, linestyle='--', color='black');\n",
    "ylabel(r\"$E_0(u)$\")\n",
    "xlabel(r\"$u=x/N$\");\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(0, -0.0014), xycoords='data',\n",
    "    xytext=(0.5/6, -0.0014), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x6 padding', xy=(0.5/12, -0.0014), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points')\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(0, -0.0012), xycoords='data',\n",
    "    xytext=(0.5/4, -0.0012), textcoords='data',\n",
    "    arrowprops={'arrowstyle': '<->'})\n",
    "plt.annotate(\n",
    "    'x4 padding', xy=(0.5/8, -0.0012), xycoords='data', ha='center',\n",
    "    xytext=(0, 1), textcoords='offset points');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 435
    },
    "id": "Nie2w5jpYzeA",
    "outputId": "a2693dc9-5203-4fe6-adc0-72d83cc1923c"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "#lab = [\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic\",\"Quintic (GB)\",\"Quintic (MZV)\",\"Quintic (JE)\",\"mks2021\"]\n",
    "lab = [\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic\",\"Quintic (GB)\",\"Quintic (MZV)\",\"Quintic (JE)\"]\n",
    "lsStyle = [\"-\"]*len(lab)\n",
    "#lsStyle[-2]=\"--\"\n",
    "lsStyle[-1]=\"-.\"\n",
    "#kerU = [hatlanczos3,hatlanczos5,hat3,hat5Gary,hat5MZV,hat5JE,hatmks2021]\n",
    "\n",
    "kerU = [hatlanczos3,hatlanczos5,hat3,hat5Gary,hat5MZV,hat5JE]\n",
    "\n",
    "\n",
    "\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "\n",
    "uv = np.linspace(0,5,10_000)\n",
    "for ik,kU in enumerate(kerU):\n",
    "  plot(uv,np.abs(kU(uv)),label=lab[ik],ls=lsStyle[ik],c=colors[ik+2])\n",
    "ylim([0,0.01])\n",
    "xlabel(r\"$u$\");\n",
    "ylabel(r\"$\\hat{K}(u)$\")\n",
    "grid();\n",
    "legend(loc='center right', bbox_to_anchor=(1.45, 0.5), fontsize=12);\n",
    "axhline(0.001, linestyle='--', color='black');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve\n",
    "\n",
    "class RootFinder:\n",
    "    def __init__(self, start, stop, step=0.01, root_dtype=\"float64\", xtol=1e-9):\n",
    "\n",
    "        self.start = start\n",
    "        self.stop = stop\n",
    "        self.step = step\n",
    "        self.xtol = xtol\n",
    "        self.roots = np.array([], dtype=root_dtype)\n",
    "\n",
    "    def add_to_roots(self, x):\n",
    "\n",
    "        if (x < self.start) or (x > self.stop):\n",
    "            return  # outside range\n",
    "        if any(abs(self.roots - x) < self.xtol):\n",
    "            return  # root already found.\n",
    "\n",
    "        self.roots = np.append(self.roots, x)\n",
    "\n",
    "    def find(self, f, *args):\n",
    "        current = self.start\n",
    "\n",
    "        for x0 in np.arange(self.start, self.stop + self.step, self.step):\n",
    "            if x0 < current:\n",
    "                continue\n",
    "            x = self.find_root(f, x0, *args)\n",
    "            if x is None:  # no root found.\n",
    "                continue\n",
    "            current = x\n",
    "            self.add_to_roots(x)\n",
    "\n",
    "        return self.roots\n",
    "\n",
    "    def find_root(self, f, x0, *args):\n",
    "\n",
    "        x, _, ier, _ = fsolve(f, x0=x0, args=args, full_output=True, xtol=self.xtol)\n",
    "        if ier == 1:\n",
    "            return x[0]\n",
    "        return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lab = [\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic    \",\"Quintic (GB)\",\"Quintic (MZV)\",\"Quintic (JE)\"]\n",
    "kerU = [hatlanczos3,hatlanczos5,hat3,hat5Gary,hat5MZV,hat5JE]\n",
    "for ik,kU in enumerate(kerU):\n",
    "    r = RootFinder(1, 5, 0.01)\n",
    "    def func(u):\n",
    "        return np.abs(kU(u)) - 0.001\n",
    "    uinit = 1.\n",
    "    umax = r.find(func)\n",
    "    print(lab[ik],'\\t',np.around(np.max(umax),2))\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lab = [\"Lanczos Norm 3\",\"Lanczos Norm 5\",\"Cubic    \",\"Quintic (GB)\",\"Quintic (MZV)\",\"Quintic (JE)\"]\n",
    "kerU = [hatlanczos3,hatlanczos5,hat3,hat5Gary,hat5MZV,hat5JE]\n",
    "for ik,kU in enumerate(kerU):\n",
    "    print(\">>>>>>> \",lab[ik],\"<<<<<<<\")\n",
    "    for thres in [1/4,1/8,1/12]:\n",
    "        uv = np.linspace(-thres,thres,50_000,endpoint=True)\n",
    "        ghost=np.max(np.abs(kU(uv+1)))\n",
    "        uv = np.linspace(0,thres,25_000,endpoint=True)\n",
    "        e0=np.max(E0(uv,kU))\n",
    "        print(np.around(thres,3),\"\\t\",np.around(e0,5),\"\\t\",np.around(ghost,5),\"\\t\",np.around(np.max((e0,ghost)),5))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare image for x->u DFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "s0-eExYsoUoh"
   },
   "source": [
    "## Make a bullseye  32x32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from google.colab import drive\n",
    "#drive.mount('/content/drive')\n",
    "#cd /content/drive/MyDrive/interpolPoly/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 452
    },
    "id": "qxNtD7oVB3kO",
    "outputId": "07bb1f7a-6fdd-44af-dbc4-2e2afe66d389"
   },
   "outputs": [],
   "source": [
    "#imgOrig = plt.imread(\"bullseye_32x32.png\")\n",
    "#imshow(imgOrig,\"gray_r\");colorbar();\n",
    "#title(\"Original image 32x32\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imgOrig = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.15294118,\n",
    "  0.40392157, 0.6431373, 0.7647059, 0.84313726, 0.84313726, 0.7647059,\n",
    "  0.6431373, 0.40392157, 0.15294118, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.003921569, 0.30588236,\n",
    "  0.7647059, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.7647059,\n",
    "  0.30588236, 0.003921569, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.15294118, 0.8117647, 1.0, 1.0,\n",
    "  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.8117647,\n",
    "  0.15294118, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.34901962, 0.9529412, 1.0, 1.0, 1.0,\n",
    "  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.9529412, 0.34901962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.4745098, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.4745098, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.30588236, 0.99607843, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  1.0, 0.78039217, 0.43529412, 0.23921569, 0.12156863, 0.12156863,\n",
    "  0.23921569, 0.43529412, 0.78039217, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.99607843, 0.30588236, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.16470589, 0.9490196, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.7882353, 0.22352941, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.22352941, 0.7882353, 1.0, 1.0, 1.0, 1.0, 1.0, 0.9490196, 0.16470589,\n",
    "  0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.7921569, 1.0, 1.0, 1.0, 1.0, 1.0, 0.627451,\n",
    "  0.02745098, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.02745098, 0.627451, 1.0, 1.0, 1.0, 1.0, 1.0, 0.7921569, 0.0, 0.0,\n",
    "  0.0,],\n",
    " [0.0, 0.0, 0.24313726, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5803922, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5803922,\n",
    "  1.0, 1.0, 1.0, 1.0, 1.0, 0.23921569, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.6901961, 1.0, 1.0, 1.0, 1.0, 0.85882354, 0.03529412, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.03529412, 0.85882354, 1.0, 1.0, 1.0, 1.0, 0.6901961, 0.0, 0.0,],\n",
    " [0.0, 0.13725491, 0.99215686, 1.0, 1.0, 1.0, 0.99607843, 0.23529412,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0.25490198, 0.7058824, 0.9098039, 0.9098039,\n",
    "  0.7058824, 0.25490198, 0.0, 0.0, 0.0, 0.0, 0.0, 0.23529412, 0.99607843,\n",
    "  1.0, 1.0, 1.0, 0.99215686, 0.13725491, 0.0,],\n",
    " [0.0, 0.38039216, 1.0, 1.0, 1.0, 1.0, 0.8117647, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.47058824, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.47058824, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.8117647, 1.0, 1.0, 1.0, 1.0, 0.38039216, 0.0,],\n",
    " [0.0, 0.5254902, 1.0, 1.0, 1.0, 1.0, 0.54901963, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.25490198, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.25490198, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.54901963, 1.0, 1.0, 1.0, 1.0, 0.5254902, 0.0,],\n",
    " [0.0, 0.6666667, 1.0, 1.0, 1.0, 1.0, 0.28627452, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.7058824, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.7058824, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.28627452, 1.0, 1.0, 1.0, 1.0, 0.6666667, 0.0,],\n",
    " [0.0, 0.8117647, 1.0, 1.0, 1.0, 1.0, 0.14901961, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.9098039, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.9098039, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.15294118, 1.0, 1.0, 1.0, 1.0, 0.80784315, 0.0,],\n",
    " [0.0, 0.80784315, 1.0, 1.0, 1.0, 1.0, 0.23529412, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.9098039, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.9098039, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.23529412, 1.0, 1.0, 1.0, 1.0, 0.80784315, 0.0,],\n",
    " [0.0, 0.6666667, 1.0, 1.0, 1.0, 1.0, 0.32156864, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.7058824, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.7058824, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.32156864, 1.0, 1.0, 1.0, 1.0, 0.6666667, 0.0,],\n",
    " [0.0, 0.52156866, 1.0, 1.0, 1.0, 1.0, 0.44705883, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.25490198, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.25490198, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.44705883, 1.0, 1.0, 1.0, 1.0, 0.52156866, 0.0,],\n",
    " [0.0, 0.38039216, 1.0, 1.0, 1.0, 1.0, 0.8392157, 0.003921569, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.47058824, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.47058824, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.003921569, 0.8392157, 1.0, 1.0, 1.0, 1.0, 0.38039216,\n",
    "  0.0,],\n",
    " [0.0, 0.13725491, 0.99215686, 1.0, 1.0, 1.0, 1.0, 0.29803923, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.25490198, 0.7058824, 0.9098039, 0.9098039, 0.7058824,\n",
    "  0.25490198, 0.0, 0.0, 0.0, 0.0, 0.0, 0.29803923, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.99215686, 0.13725491, 0.0,],\n",
    " [0.0, 0.0, 0.6901961, 1.0, 1.0, 1.0, 1.0, 0.78039217, 0.019607844, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.019607844, 0.78431374, 1.0, 1.0, 1.0, 1.0, 0.6901961, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.24313726, 1.0, 1.0, 1.0, 1.0, 1.0, 0.64705884, 0.003921569,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.003921569, 0.64705884, 1.0, 1.0, 1.0, 1.0, 1.0, 0.23921569, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.7921569, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5803922,\n",
    "  0.015686275, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.015686275, 0.5803922, 1.0, 1.0, 1.0, 1.0, 1.0, 0.7921569, 0.0, 0.0,\n",
    "  0.0,],\n",
    " [0.0, 0.0, 0.0, 0.16470589, 0.9490196, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.81960785, 0.23137255, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.23137255, 0.81960785, 1.0, 1.0, 1.0, 1.0, 1.0, 0.9490196, 0.16470589,\n",
    "  0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.30588236, 0.99607843, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.99607843, 0.76862746, 0.45490196, 0.22352941, 0.1254902, 0.1254902,\n",
    "  0.22352941, 0.45490196, 0.76862746, 0.99607843, 1.0, 1.0, 1.0, 1.0,\n",
    "  1.0, 0.99607843, 0.30588236, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.4745098, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.4745098, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.34901962, 0.9529412, 1.0, 1.0, 1.0,\n",
    "  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "  0.9529412, 0.34901962, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.15294118, 0.8117647, 1.0, 1.0,\n",
    "  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.8117647,\n",
    "  0.15294118, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.003921569, 0.30588236,\n",
    "  0.7647059, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.7647059,\n",
    "  0.30588236, 0.003921569, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.15294118,\n",
    "  0.40392157, 0.6431373, 0.7647059, 0.84313726, 0.84313726, 0.7647059,\n",
    "  0.6431373, 0.40392157, 0.15294118, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0,],\n",
    " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "  0.0, 0.0, 0.0, 0.0,]],dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(imgOrig,\"gray_r\");colorbar();\n",
    "title(\"Original image 32x32\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "id": "EUdKMeI_VmwG",
    "outputId": "6f66fa3e-e42b-496f-a2ba-408e93637bc2"
   },
   "outputs": [],
   "source": [
    "#scatter(np.arange(32),imgOrig[:,16]);\n",
    "#plot(imgOrig[:,16]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "B5l9L6TBj0Nz"
   },
   "outputs": [],
   "source": [
    "N1o,N2o=imgOrig.shape\n",
    "N1o,N2o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xkzcIDKk8f_X"
   },
   "source": [
    "## (NOT USED) From the original 32x32 image, makes interpolation on a x2 oversampled grid (lancsoz-3 normalized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "LiScNR3-LELM"
   },
   "outputs": [],
   "source": [
    "def interpX(img,x1,x2,ker):\n",
    "  N1,N2 = img.shape\n",
    "  res=[ker(x-p) for p in range(-N1//2,N1//2,1) for x in x1]\n",
    "  mtxKx = np.array(res).reshape((N1,-1))\n",
    "  res=[ker(x-p) for p in range(-N2//2,N2//2,1) for x in x2]\n",
    "  mtxKy = np.array(res).reshape((N2,-1))\n",
    "\n",
    "  # mtxKx.T .       I(N1,N2) . mtxKy\n",
    "  #print(mtxKx.T.shape,img.shape,mtxKy.shape)\n",
    "  return mtxKx.T @ img @  mtxKy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "9lsWWS5Xsi2d",
    "outputId": "16a0d688-7b5b-459f-e119-20de413b1a5c"
   },
   "outputs": [],
   "source": [
    "N1ov = 2*N1o\n",
    "N2ov = 2*N2o\n",
    "N1o,N2o,N1ov,N2ov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "KDOnOUJbrD-7"
   },
   "outputs": [],
   "source": [
    "#x1 = np.arange(0,N1o,N1o/N1ov)#   np.linspace(0,N1,2*N1,endpoint=False)   # care with linespace\n",
    "#x2 = np.arange(0,N2o,N2o/N2ov)#   np.linspace(0,N2,2*N2,endpoint=False)\n",
    "#test = interp(img,x1,x2,h3)\n",
    "#test.min(), test.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 452
    },
    "id": "8KMgU8n4Ll3r",
    "outputId": "d1c8a6bb-54f5-4240-e8ef-ae72f143c5fc"
   },
   "outputs": [],
   "source": [
    "x1 = np.arange(-N1o//2,N1o//2,N1o/N1ov)\n",
    "x2 = np.arange(-N2o//2,N2o//2,N1o/N1ov)\n",
    "orig_x2 = interpX(imgOrig,x1,x2,lanczos3)  \n",
    "imshow(orig_x2,\"gray_r\"); colorbar();\n",
    "title(\"Original image 2x oversampled (Lanczos norm 3rd)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "id": "WF02a2KFW5bY",
    "outputId": "edcfe237-ef69-4a39-eb6d-86498f2a0e3f"
   },
   "outputs": [],
   "source": [
    "scatter(np.arange(64),orig_x2[:,32]);\n",
    "plot(orig_x2[:,32]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DHozknNaoN_V"
   },
   "source": [
    "##  Pad this 32x32 *original image* (not the x2 interpolated) to get a 128x128 image before x->u DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ZGAYUv4Qg5PD"
   },
   "outputs": [],
   "source": [
    "orig_pad_128 = np.pad(imgOrig, (128-32)//2, pad_with, padder=0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Jz870rjchOtj",
    "outputId": "8e88870c-767d-44f5-8700-05b9a281b284"
   },
   "outputs": [],
   "source": [
    "orig_pad_128.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DIaxaHQ5hFtC"
   },
   "outputs": [],
   "source": [
    "imshow(orig_pad_128,\"gray_r\"); colorbar();\n",
    "title(\"Original image padded to 128x128\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "id": "WfKK5lJxYKGk",
    "outputId": "06e80758-91f8-4373-ba33-d626eb934dd2"
   },
   "outputs": [],
   "source": [
    "#scatter(np.arange(128),orig_pad_128[:,64]);\n",
    "#plot(orig_pad_128[:,64]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "37g9mY5mAomU"
   },
   "source": [
    "# Find DFT in 1D & 2D as defined by Gary et al. working with np.fft framework"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J8ZC8Em-0gQ2"
   },
   "outputs": [],
   "source": [
    "x=np.random.normal(size=(10,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "GvU38-VG5FIV"
   },
   "outputs": [],
   "source": [
    "def dft1D(x):\n",
    "  # mimic the numpy FFT\n",
    "  N = x.shape[0]\n",
    "  mtx = np.array([np.exp(-2.0*1j*np.pi * m * k/N) for m in range(N) for k in range(N)]).reshape(N,-1)\n",
    "  #print(mtx.shape, N)\n",
    "  return mtx @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "juhThZuzIH8m"
   },
   "outputs": [],
   "source": [
    "def invdft1D(x):\n",
    "  # mimic the numpy iFFT\n",
    "  N = x.shape[0]\n",
    "  mtx = np.array([np.exp(2.0*1j*np.pi * m * k/N) for m in range(N) for k in range(N)]).reshape(N,-1)\n",
    "  return (mtx @ x)/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cyyvpFR4_P3g"
   },
   "outputs": [],
   "source": [
    "def dft1Dshift(x):\n",
    "  # mimic the Gary et al FFT\n",
    "  N = x.shape[0]\n",
    "  mtx = np.array([np.exp(-2.0*1j*np.pi * m * k/N) for m in range(-N//2,N//2) for k in range(-N//2,N//2)]).reshape(N,-1)\n",
    "  #print(mtx.shape, N)\n",
    "  return mtx @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5GAFdGJeMOWU"
   },
   "outputs": [],
   "source": [
    "def npdft1Dshift(x):\n",
    "  # mimic Gary et al 1D DFT using numpy FFT\n",
    "  return np.fft.fftshift(np.fft.fft(x))*np.array([(-1)**k for k in range(-x.shape[0]//2,x.shape[0]//2)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nGmOcEBsIbri"
   },
   "outputs": [],
   "source": [
    "def invdft1Dshift(x):\n",
    "  # mimic the Gary et al iFFT\n",
    "  N = x.shape[0]\n",
    "  mtx = np.array([np.exp(2.0*1j*np.pi * m * k/N) for m in range(-N//2,N//2) for k in range(-N//2,N//2)]).reshape(N,-1)\n",
    "  #print(mtx.shape, N)\n",
    "  return (mtx @ x)/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bT-MfOa-OgVQ"
   },
   "outputs": [],
   "source": [
    "def npinvdft1Dshift(x):\n",
    "  # mimic Gary et al 1D i DFT using numpy iFFT\n",
    "  return np.fft.ifftshift(np.fft.ifft(x))*np.array([(-1)**k for k in range(-x.shape[0]//2,x.shape[0]//2)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "QL1w7BPPH5dE",
    "outputId": "deda8122-6b67-4f4d-a784-ef37b7f68535"
   },
   "outputs": [],
   "source": [
    "np.allclose(dft1D(x),np.fft.fft(x))  # peddestrian DFT mimic numpy DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "ZfGmvX2RPx8A",
    "outputId": "d6dd2669-7630-4b07-c7bf-3552cd0e05ee"
   },
   "outputs": [],
   "source": [
    "np.allclose(invdft1D(x),np.fft.ifft(x)) # peddestrian iDFT mimic numpy iDFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "tHoy3J3vMzS8",
    "outputId": "0f9d45ae-72d2-410d-db4c-db3833fb26c1"
   },
   "outputs": [],
   "source": [
    "np.allclose(dft1Dshift(x),npdft1Dshift(x)) # numpy FFT mimic Gary DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "2fLhYJNLPUj1",
    "outputId": "2f17cf50-daa2-4fc5-c5b9-da01d43385b5"
   },
   "outputs": [],
   "source": [
    "np.allclose(npinvdft1Dshift(x),invdft1Dshift(x))# numpy iFFT mimic Gary iDFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "K9OSvmfxCCXd",
    "outputId": "cfd6e958-2e12-4e8d-920a-ccfbfc0e8ed5"
   },
   "outputs": [],
   "source": [
    "np.allclose(invdft1D(dft1D(x)),x) # self consistance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "T-zNZGYTIpXS",
    "outputId": "a727bb43-8cc3-460a-8453-b52154de6d36"
   },
   "outputs": [],
   "source": [
    "np.allclose(invdft1Dshift(dft1Dshift(x)),x) # self consistance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "eGNc5ZOOODoS",
    "outputId": "c229e917-4236-4066-d01c-3e5b1c9d08f6"
   },
   "outputs": [],
   "source": [
    "np.allclose(npinvdft1Dshift(npdft1Dshift(x)),x) # self consistance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5xtQERvOOxOg"
   },
   "outputs": [],
   "source": [
    "x=np.random.normal(size=(4,6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Gt8pgJe-BRfy",
    "outputId": "215bd004-3f95-4b7d-c7f0-e5dac72e3a47"
   },
   "outputs": [],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "PPNvkTzyBDFs"
   },
   "outputs": [],
   "source": [
    "def dft2D(x):\n",
    "  # mimic the numpy FFT\n",
    "  N1,N2 = x.shape\n",
    "  mtx1 = np.array([np.exp(-2.0*1j*np.pi * n1 * k1/N1) for n1 in range(N1) for k1 in range(N1)]).reshape(N1,-1)\n",
    "  mtx2 = np.array([np.exp(-2.0*1j*np.pi * n2 * k2/N2) for n2 in range(N2) for k2 in range(N2)]).reshape(N2,-1)\n",
    "  return  mtx1.T @ x @ mtx2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "L3TqVmIeMNNE",
    "outputId": "33c41062-d53d-446f-cc13-c517cda8d348"
   },
   "outputs": [],
   "source": [
    "np.allclose(dft2D(x),np.fft.fft2(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cVQZKxdzMwND"
   },
   "outputs": [],
   "source": [
    "def dft2Dshift(x):\n",
    "  # mimic the Gary et al FFT\n",
    "  N1,N2 = x.shape\n",
    "  mtx1 = np.array([np.exp(-2.0*1j*np.pi * n1 * k1/N1) for n1 in range(-N1//2,N1//2) for k1 in range(-N1//2,N1//2)]).reshape(N1,-1)\n",
    "  mtx2 = np.array([np.exp(-2.0*1j*np.pi * n2 * k2/N2) for n2 in range(-N2//2,N2//2) for k2 in range(-N2//2,N2//2)]).reshape(N2,-1)\n",
    "  return  mtx1.T @ x @ mtx2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dWilCayGMyax"
   },
   "outputs": [],
   "source": [
    "def npdft2Dshift(x):\n",
    "  # mimic Gary et al 2D DFT using numpy FFT\n",
    "  N1,N2 = x.shape\n",
    "  mtxsig = np.array([(-1)**(k1+k2) for k1 in range(-N1//2,N1//2) for k2 in range(-N2//2,N2//2)]).reshape(N1,N2)\n",
    "  return np.fft.fftshift(np.fft.fft2(x)) * mtxsig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "B3-BtzksPbvc",
    "outputId": "c261b383-153a-4b8c-b0c4-966263c46d32"
   },
   "outputs": [],
   "source": [
    "np.allclose(dft2Dshift(x),npdft2Dshift(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fmIXIBukPi7X"
   },
   "outputs": [],
   "source": [
    "def invdft2D(x):\n",
    "  # mimic the numpy iFFT\n",
    "  N1,N2 = x.shape\n",
    "  mtx1 = np.array([np.exp(2.0*1j*np.pi * n1 * k1/N1) for n1 in range(N1) for k1 in range(N1)]).reshape(N1,-1)\n",
    "  mtx2 = np.array([np.exp(2.0*1j*np.pi * n2 * k2/N2) for n2 in range(N2) for k2 in range(N2)]).reshape(N2,-1)\n",
    "  return  (mtx1.T @ x @ mtx2)/(N1*N2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "NvdObBCsRAFe",
    "outputId": "f7107c42-640d-4c96-c97d-ad3f204284e9"
   },
   "outputs": [],
   "source": [
    "np.allclose(invdft2D(x),np.fft.ifft2(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "pCAO8tXZRFUE"
   },
   "outputs": [],
   "source": [
    "def invdft2Dshift(x):\n",
    "  # mimic the Gary et al iFFT\n",
    "  N1,N2 = x.shape\n",
    "  mtx1 = np.array([np.exp(2.0*1j*np.pi * n1 * k1/N1) for n1 in range(-N1//2,N1//2) for k1 in range(-N1//2,N1//2)]).reshape(N1,-1)\n",
    "  mtx2 = np.array([np.exp(2.0*1j*np.pi * n2 * k2/N2) for n2 in range(-N2//2,N2//2) for k2 in range(-N2//2,N2//2)]).reshape(N2,-1)\n",
    "  return  (mtx1.T @ x @ mtx2)/(N1*N2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "_jqgmYIWRbZW",
    "outputId": "6f504295-ac50-4980-d9ee-e0c62f9e8246"
   },
   "outputs": [],
   "source": [
    "np.allclose(invdft2Dshift(dft2Dshift(x)),x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "toT_eqGwRk6I"
   },
   "outputs": [],
   "source": [
    "def npinvdft2Dshift(x):\n",
    "  # mimic Gary et al 2D DFT using numpy FFT\n",
    "  N1,N2 = x.shape\n",
    "  mtxsig = np.array([(-1)**(k1+k2) for k1 in range(-N1//2,N1//2) for k2 in range(-N2//2,N2//2)]).reshape(N1,N2)\n",
    "  return np.fft.ifftshift(np.fft.ifft2(x)) * mtxsig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "6JU2s4bARwXr",
    "outputId": "a037c488-9729-4c63-be1a-bdf90eb0c274"
   },
   "outputs": [],
   "source": [
    "np.allclose(invdft2Dshift(x),npinvdft2Dshift(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "bH9HWClmXVXB",
    "outputId": "4deca66b-4e1c-4b34-aca8-2e7f3b1852b1"
   },
   "outputs": [],
   "source": [
    "np.allclose(npinvdft2Dshift(npdft2Dshift(x)),x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EfZSRUlKSBqU",
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "So, we can use npdft2Dshift and npinvdft2Dshift to perform Gary et al 2D FFT and inverse 2D FFT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eJMZM_V4TXfB"
   },
   "source": [
    "# Fourier space Interpolation using wrapped kernels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrapped kernels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ito4mkmXR6zz"
   },
   "outputs": [],
   "source": [
    "def KuSincWrapped(u,N):\n",
    "  # \"sinc kernel wrapped\"\n",
    "  return np.exp(1j * np.pi * u) * np.sinc(N*u)/np.sinc(u) # formule (10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 474
    },
    "id": "3LEzZmJPPefr",
    "outputId": "6aaea1db-750a-4e85-f1a7-2a8c3521f36b"
   },
   "outputs": [],
   "source": [
    "u=np.linspace(-2.5,2.5,10_000)\n",
    "for n in [5,10,100]:\n",
    "  plot(u,np.abs(np.sinc(n*u)/np.sinc(u)),label=f\"N={n}\")\n",
    "xlabel(\"x\")\n",
    "legend()\n",
    "title(r\"$|sinc(Nx)/sinc(x)|$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lanczos_p(x,m,N):\n",
    "  x = np.abs(x)\n",
    "  th=m/N\n",
    "  return np.piecewise(x, [x<th,x>=th],[lambda x: lanczos(N*x,m)/lanczos(x,m), lambda x:0])\n",
    "\n",
    "def rlanczosWrapped(u,N,m=3):\n",
    "    return np.sign(0.5-np.floor(u+0.5)%2) * lanczos_p(u-np.floor(u+0.5),m,N)\n",
    "\n",
    "def lanczosWrapped(u,N,m=3):\n",
    "    return np.exp(1j * np.pi * u) * rlanczosWrapped(u,N,m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cubic_p(x,N):\n",
    "  x = np.abs(x)\n",
    "  th=2./N\n",
    "  return np.piecewise(x, [x<th,x>=th],[lambda x: h3(N*x)/h3(x), lambda x:0])\n",
    "def rcubicWrapped(u,N):\n",
    "    return np.sign(0.5-np.floor(u+0.5)%2)*cubic_p(u-np.floor(u+0.5),N)\n",
    "def cubicWrapped(u,N):\n",
    "    return np.exp(1j * np.pi * u) * rcubicWrapped(u,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quinticBG_p(x,N):\n",
    "  x = np.abs(x)\n",
    "  th=3./N\n",
    "  return np.piecewise(x, [x<th,x>=th],[lambda x: h5Gary(N*x)/h5Gary(x), lambda x:0])\n",
    "def rquinticBGWrapped(u,N):\n",
    "    return np.sign(0.5-np.floor(u+0.5)%2)*quinticBG_p(u-np.floor(u+0.5),N)\n",
    "def quinticBGWrapped(u,N):\n",
    "    return np.exp(1j * np.pi * u) * rquinticBGWrapped(u,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quinticJE_p(x,N):\n",
    "  x = np.abs(x)\n",
    "  th=3./N\n",
    "  return np.piecewise(x, [x<th,x>=th],[lambda x: h5JE(N*x)/h5JE(x), lambda x:0])\n",
    "def rquinticJEWrapped(u,N):\n",
    "    return np.sign(0.5-np.floor(u+0.5)%2)*quinticJE_p(u-np.floor(u+0.5),N)\n",
    "def quinticJEWrapped(u,N):\n",
    "    return np.exp(1j * np.pi * u) * rquinticJEWrapped(u,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "u=np.linspace(-2.5,2.5,10_000)\n",
    "N = 32\n",
    "plot(u,np.sinc(N*u)/np.sinc(u),label=\"Periodic Sinc\")\n",
    "plot(u,rlanczosWrapped(u,N,m=3),label=\"Periodic Lanczos 3\",ls='--')\n",
    "plot(u,rlanczosWrapped(u,N,m=5),label=\"Periodic Lanczos 5\",ls='--')\n",
    "plot(u,rcubicWrapped(u,N),label=\"Periodic Cubic\",ls='--')\n",
    "plot(u,rquinticBGWrapped(u,N),label=\"Periodic Quintic BG\",ls='--')\n",
    "plot(u,rquinticJEWrapped(u,N),label=\"Periodic Quintic JE\",ls='--')\n",
    "legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "u=np.linspace(-0.2,0.2,10_000)\n",
    "N = 32\n",
    "plot(u,np.sinc(N*u)/np.sinc(u),label=\"Periodic Sinc\")\n",
    "plot(u,rlanczosWrapped(u,N,m=3),label=\"Periodic Lanczos 3\",ls='--')\n",
    "plot(u,rlanczosWrapped(u,N,m=5),label=\"Periodic Lanczos 5\",ls='--')\n",
    "plot(u,rcubicWrapped(u,N),label=\"Periodic Cubic\",ls='--')\n",
    "plot(u,rquinticBGWrapped(u,N),label=\"Periodic Quintic BG\",ls='--')\n",
    "plot(u,rquinticJEWrapped(u,N),label=\"Periodic Quintic JE\",ls='--')\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paper Eq 9 in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NaIEwfZXzPlC"
   },
   "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,
   "metadata": {
    "id": "VpzeUT1xncqT"
   },
   "outputs": [],
   "source": [
    "def fhat2DLanczos3(img,u1,u2):\n",
    "  raise \"Not to be used\" #use the generic fhatKernel\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([lanczosWrapped(u-k1/N1,N1) for k1 in range(-N1//2,N1//2) for u in u1]).reshape(N1,-1)\n",
    "  #mtx2 = np.array([lanczosWrapped(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,
   "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": "code",
   "execution_count": null,
   "metadata": {
    "id": "ibEeMP_94dLO"
   },
   "outputs": [],
   "source": [
    "def fhat2D(img,u1,u2,ker):\n",
    "  raise \"Not to be used\"\n",
    "  # apply DFT to img\n",
    "  #N1,N2 = img.shape\n",
    "  #def cplx_ker(v,N):\n",
    "  #  sum = ker(N*v)\n",
    "  #  sign = -1\n",
    "  #  for i in range(1,2): # Todo better\n",
    "  #    sum +=  sign*(ker(N*(v+i))+ker(N*(v-i)))\n",
    "  #    sign *= -1\n",
    "  #  return np.exp(1j * np.pi * v)*sum\n",
    "  #atildekl = npdft2Dshift(img) # formule (6) of fourier coeff in 2D\n",
    "  #mtx1 = np.array([cplx_ker(u-k1/N1,N1) for k1 in range(-N1//2,N1//2) for u in u1]).reshape(N1,-1)\n",
    "  #mtx2 = np.array([cplx_ker(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",
   "metadata": {},
   "source": [
    "## Generic TF[Kx] kernel matrix in 2D Fourier grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "h2DCeOXVbh04"
   },
   "outputs": [],
   "source": [
    "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,
   "metadata": {
    "id": "dXlUyJ9Iy8Pf"
   },
   "outputs": [],
   "source": [
    "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the new u-grid for Fourier resampling "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "OL1pUuSHIvYS",
    "outputId": "e79f3ad3-4d7a-47e9-9378-a50e71c982cf"
   },
   "outputs": [],
   "source": [
    "N1pad,N2pad = orig_pad_128.shape\n",
    "N1pad,N2pad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "_kejIJiai6Mz",
    "outputId": "75a4a6a7-df9f-4bab-eac7-39279dfb2252"
   },
   "outputs": [],
   "source": [
    "# hypo: N1pad=N2pad\n",
    "M=N1pad*4 #  width/height of the final image (Figure 2 if bulleteyes is 32x32)\n",
    "dx=1\n",
    "du = 1/(M*dx)\n",
    "umax = (M/2)*du\n",
    "umax,du"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6_e0BB7Z3Jce"
   },
   "outputs": [],
   "source": [
    "u1=np.arange(-umax,umax,du)\n",
    "u2=np.arange(-umax,umax,du)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "nlza6ubhxuHV",
    "outputId": "df4773e5-f62f-4737-ac64-605bdfe96e28"
   },
   "outputs": [],
   "source": [
    "u1.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "use Lanczos norm 3rd order for Kx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Kxhat_mtx = Kxhat2D(u1,u2,hatlanczos3)  # Lanczos norm 3rd order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ku is the Sinc-periodic kernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TF[fd] Eq 7-9 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hGIW_Qgu-I8E"
   },
   "outputs": [],
   "source": [
    "fhat_mtx_ideal = fhat2DSincWrapped(orig_pad_128,u1,u2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "x8bhHH2q-RL6"
   },
   "outputs": [],
   "source": [
    "imshow(np.abs(fhat_mtx_ideal));colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TF\\[F\\](u) Eq. 11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qdR_rnBMGy-v"
   },
   "outputs": [],
   "source": [
    "Fhat_mtx_ideal = Kxhat_mtx * fhat_mtx_ideal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### F[x] as inverse transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jkLKZlJIHDjR"
   },
   "outputs": [],
   "source": [
    "F_real_ideal = npinvdft2Dshift(Fhat_mtx_ideal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### cleaning small numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "FUQg8JI7HLxt"
   },
   "outputs": [],
   "source": [
    "F_real_ideal = clean_tiny_cplx(F_real_ideal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "BQsYI38BHUf9"
   },
   "outputs": [],
   "source": [
    "F_real_ideal[np.abs(F_real_ideal)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0OodKkA1ucAh"
   },
   "outputs": [],
   "source": [
    "imshow(F_real_ideal,\"gray_r\",vmin=-0.003,vmax=0.003);colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ku = lanczos3\n",
    "Use the same kind of procedure as for Sinc-periodic but with lanczos3 as Ku. (Eq. 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "envB_-wRZnrG"
   },
   "outputs": [],
   "source": [
    "#fhat_mtx = fhat2D(orig_pad_128,u1,u2,lanczos3)  # Ku = lanczos3\n",
    "fhat_lanczos3_mtx = fhatKernel(orig_pad_128,u1,u2,lanczosWrapped) # Ku = lanczos3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "xLaaUv8O2tWP"
   },
   "outputs": [],
   "source": [
    "#imshow(np.abs(fhat_lanczos3_mtx));colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Vpkotu_F7mdE"
   },
   "outputs": [],
   "source": [
    "f_lanczos3_real = npinvdft2Dshift(fhat_lanczos3_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "UZSvg1m87tRo"
   },
   "outputs": [],
   "source": [
    "f_lanczos3_real =clean_tiny_cplx(f_lanczos3_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jNJ4GOni76YY"
   },
   "outputs": [],
   "source": [
    "imshow(f_lanczos3_real,\"gray_r\");colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lsCaca5CfM8F"
   },
   "outputs": [],
   "source": [
    "plot(np.abs(f_lanczos3_real[M//2,:]))\n",
    "yscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fI3EUbKe24PG"
   },
   "outputs": [],
   "source": [
    "Fhat_lanczos3_mtx = Kxhat_mtx * fhat_lanczos3_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "sOv26R9u5p0b"
   },
   "outputs": [],
   "source": [
    "#imshow(np.abs(Fhat_mtx));colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "x4eHeVVG5rx2"
   },
   "outputs": [],
   "source": [
    "F_lanczos3_real = npinvdft2Dshift(Fhat_lanczos3_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "XI_JqDAa6lAk"
   },
   "outputs": [],
   "source": [
    "F_lanczos3_real = clean_tiny_cplx(F_lanczos3_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "opCQ8oUg6v6Z"
   },
   "outputs": [],
   "source": [
    "F_lanczos3_real[np.abs(F_lanczos3_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WqEyXvC_684j"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(5,5))\n",
    "imshow(F_lanczos3_real,\"gray_r\",vmin=-0.003,vmax=0.003);colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "MRyxONIFtEA7"
   },
   "outputs": [],
   "source": [
    "figure(figsize=(5,5))\n",
    "imshow(np.log10(np.abs(F_lanczos3_real)+1e-12),\"jet\",vmin=-6);colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "i9t2b044t4al"
   },
   "outputs": [],
   "source": [
    "imshow(crop_at(F_lanczos3_real,64,64),\"gray_r\");colorbar();\n",
    "title(\"Central image of F(x) (Lanzcos 3 in real & Fourier space)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "noLV9VCAsZ9-"
   },
   "outputs": [],
   "source": [
    "imshow(crop_at(F_lanczos3_real,64,64,128,0),\"gray_r\");colorbar();\n",
    "title(\"First ghost (128,0) of F(x) (Lanzcos 3 in real & Fourier space)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9Pi7xXSjIHNA"
   },
   "source": [
    "## Ku = Cubic\n",
    "replace lancsoz-3 Ku by Cubic kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_cubic_mtx = fhatKernel(orig_pad_128,u1,u2,cubicWrapped) # Ku = cubic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#imshow(np.abs(fhat_cubic_mtx));colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_cubic_real = npinvdft2Dshift(fhat_cubic_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_cubic_real =clean_tiny_cplx(f_cubic_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(np.abs(f_cubic_real[M//2,:]))\n",
    "yscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_cubic_mtx = Kxhat_mtx * fhat_cubic_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_cubic_real = npinvdft2Dshift(Fhat_cubic_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_cubic_real = clean_tiny_cplx(F_cubic_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_cubic_real[np.abs(F_cubic_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure(figsize=(5,5))\n",
    "imshow(F_cubic_real,\"gray_r\",vmin=-0.003,vmax=0.003);colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_cubic_real,64,64,128,0),\"gray_r\");colorbar();\n",
    "title(\"First ghost (128,0) of F(x) (Lanzcos 3 in real & Cubic Fourier space)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ku = Quintic (Gary)\n",
    "replace lancsoz-3 Ku by Gary's Quintic  kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_quinticBG_mtx = fhatKernel(orig_pad_128,u1,u2,quinticBGWrapped) # Ku = Quintic BG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_quinticBG_real = npinvdft2Dshift(fhat_quinticBG_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_quinticBG_real =clean_tiny_cplx(f_quinticBG_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(np.abs(f_quinticBG_real[M//2,:]))\n",
    "yscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_quinticBG_mtx = Kxhat_mtx * fhat_quinticBG_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticBG_real = npinvdft2Dshift(Fhat_quinticBG_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticBG_real = clean_tiny_cplx(F_quinticBG_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticBG_real[np.abs(F_quinticBG_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure(figsize=(5,5))\n",
    "imshow(F_quinticBG_real,\"gray_r\",vmin=-0.003,vmax=0.003);colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_quinticBG_real,64,64,128,0),\"gray_r\");colorbar();\n",
    "title(\"First ghost (128,0) of F(x) (Lanzcos 3 in real & Quintic BG Fourier space)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ku = Quintic (JE)\n",
    "replace lancsoz-3 Ku by mine Quintic  kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fhat_quinticJE_mtx = fhatKernel(orig_pad_128,u1,u2,quinticJEWrapped) # Ku = Quintic JE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_quinticJE_real = npinvdft2Dshift(fhat_quinticJE_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(np.abs(f_quinticJE_real[M//2,:]))\n",
    "yscale(\"log\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fhat_quinticJE_mtx = Kxhat_mtx * fhat_quinticJE_mtx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticJE_real = npinvdft2Dshift(Fhat_quinticJE_mtx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticJE_real = clean_tiny_cplx(F_quinticJE_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_quinticJE_real[np.abs(F_quinticJE_real)<1e-12]=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure(figsize=(5,5))\n",
    "imshow(F_quinticJE_real,\"gray_r\",vmin=-0.003,vmax=0.003);colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_quinticJE_real,64,64,128,0),\"gray_r\");colorbar();\n",
    "title(\"First ghost (128,0) of F(x) (Lanzcos 3 in real & Quintic JE Fourier space)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "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,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmin = -0.003 #-0.003\n",
    "vmax = 0.003  #0.003"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(F_lanczos3_real-F_real_ideal,\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Lanczos-3 norm - Ideal\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_lanczos3_real-F_real_ideal,64,64),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Lanczos-3 norm - Ideal central\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_cubic_real-F_real_ideal,64,64),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Cubic - Ideal central\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_quinticBG_real-F_real_ideal,64,64),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Quintic BG - Ideal central\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_quinticJE_real-F_real_ideal,64,64),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Quintic JE - Ideal central\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_lanczos3_real-F_real_ideal,64,64,128,0),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Lanczos-3 norm - Ideal 1st ghost\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_cubic_real-F_real_ideal,64,64,128,0),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Cubic - Ideal 1st ghost\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_quinticBG_real-F_real_ideal,64,64,128,0),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Quintic BG - Ideal 1st ghost\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(crop_at(F_quinticJE_real-F_real_ideal,64,64,128,0),\"gray_r\",vmin=vmin, vmax=vmax);colorbar();\n",
    "title(\"Quintic JE - Ideal 1st ghost\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sand Box"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyMgQcjCuWbtDGtCrHVTE61v",
   "include_colab_link": true,
   "provenance": [],
   "toc_visible": true
  },
  "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": 4
}