{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# SAR image interpolation\n",
    "\n",
    "A backprojected pseudo-polar SAR image is not an absolute picture of the ground it is referenced to an origin, the mean antenna phase centre (APC) of the platform during the aperture. The range $r$ and angle $\\theta$ of every pixel, and the carrier phase $\\exp{\\left(-j 4\\pi f_c r / c\\right)}$ stored in it measured from the origin.\n",
    "\n",
    "Two acquisitions flown from slightly different positions therefore land on two different polar grids. Before they can be compared, for example interferometry, they must be brought onto a common grid centred on a single origin. The brute force way is to backproject each pass again from its raw data onto the shared grid, but that is expensive.\n",
    "\n",
    "`torchbp.ops.polar_interp` does the same job directly in the image domain. It resamples an already-formed polar image from its old origin/grid onto a new origin/grid, correcting both:\n",
    "\n",
    "  * the geometry, each pixel's range and angle change when viewed from the new origin\n",
    "  * the carrier phase, the phase must be updated for the new range\n",
    "\n",
    "This is the core building block of repeat-pass coregistration and of fast-factorized backprojection (FFBP), where small sub-aperture images are recuresively merged into larger more detailed image.\n",
    "\n",
    "This example demonstrates `polar_interp` with synthetic data. We take one pass of raw data and form it into a polar image two ways: Directly at the new position, and by interpolating from nearby position. If `polar_interp` is correct the two results agree in amplitude and in phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "import torchbp\n",
    "import scipy.signal as signal\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "Use CUDA if available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "if torch.cuda.is_available():\n",
    "    device = \"cuda\"\n",
    "else:\n",
    "    device = \"cpu\"\n",
    "print(\"Device:\", device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Radar and acquisition geometry\n",
    "\n",
    "A small FMCW radar flown along the $y$ axis, looking towards $+x$. We simulate two passes:\n",
    "\n",
    "  * master (pass 0), the reference geometry\n",
    "  * slave (pass 1), flown from a slightly different position. The slave is offset `cross_track` metres in $x$, `along_track` metres in $y$, and `dh` metres higher.\n",
    "\n",
    "The vector between the two antenna phase centres is the baseline. Re-centring the slave image onto the master grid is exactly an origin shift by that baseline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 5.9e9          # RF center frequency (Hz)\n",
    "bw = 300e6          # Sweep bandwidth (Hz)\n",
    "tsweep = 200e-6     # Sweep length (s)\n",
    "fs = 10e6           # Sampling frequency (Hz)\n",
    "nsamples = int(fs * tsweep)   # Time-domain samples per sweep\n",
    "nsweeps = 384       # Number of azimuth measurements (array elements)\n",
    "wl = 3e8 / fc       # Wavelength\n",
    "\n",
    "altitude = 50.0     # Master platform altitude (m)\n",
    "cross_track = 1.0   # Slave offset across track, x (m)\n",
    "along_track = 0.5   # Slave offset along track, y (m)\n",
    "dh = 0.25            # Slave height above master (m)\n",
    "element_spacing = 0.25 * wl   # Azimuth sampling along the track"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## Polar imaging grid\n",
    "\n",
    "`theta` is the sine of the azimuth angle. The number of range and azimuth bins is set from the range resolution and the synthetic aperture so the grid is slightly oversampled. `signal_oversample` is how many grid range bins span one resolution cell, it appears again below when we centre the range spectrum for interpolation.\n",
    "\n",
    "Higher grid oversampling factor reduces the interpolation error. If the image is aliased (oversample < 1), then interpolation can't give the correct result as the aliased and desired sprectum overlap. This is important to get correct for interpolation to be possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "r0, r1 = 100.0, 200.0      # Ground-range extent (m)\n",
    "theta_limit = 0.7         # Azimuth half-extent (sine of angle)\n",
    "\n",
    "res = 3e8 / (2 * bw)      # Slant-range resolution\n",
    "res_interp = 2         # Grid oversampling factor\n",
    "nr = int(res_interp * (r1 - r0) / res)\n",
    "ntheta = int(1 + nsweeps * res_interp * (element_spacing / wl) * theta_limit / 0.25)\n",
    "\n",
    "grid_polar = {\"r\": (r0, r1), \"theta\": (-theta_limit, theta_limit),\n",
    "              \"nr\": nr, \"ntheta\": ntheta}\n",
    "\n",
    "dr_grid = (r1 - r0) / nr\n",
    "signal_oversample = res / dr_grid # How many pixels in image per resolution cell in image\n",
    "print(f\"nr={nr}, ntheta={ntheta}, signal_oversample={signal_oversample:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Scene\n",
    "\n",
    "A field of random complex scatterers (speckle) on a flat ground plane. We use dense grid of scatterers so that there is at least one scatterer per resolution cell, for nice coherence plot later. This is also very difficult scene for interpolation, since the whole image is filled with scatterers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.manual_seed(42)\n",
    "\n",
    "grid_proj = {\"x\": (30, 210), \"y\": (-210, 210), \"nx\": 384, \"ny\": 384}\n",
    "img = torch.randn([grid_proj[\"nx\"], grid_proj[\"ny\"]], dtype=torch.complex64, device=device)\n",
    "dem = torch.zeros([grid_proj[\"nx\"], grid_proj[\"ny\"]], dtype=torch.float32, device=device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Platform trajectories\n",
    "\n",
    "Both passes fly straight along $y$. The slave is rigidly shifted by the baseline. We weight the aperture with a Taylor window before computing each pass's phase centre, matching the window applied during image formation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_coords = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps) * element_spacing\n",
    "\n",
    "# Master (pass 0)\n",
    "pos0 = torch.zeros([nsweeps, 3], dtype=torch.float32, device=device)\n",
    "pos0[:, 1] = y_coords.to(device)\n",
    "pos0[:, 2] = altitude\n",
    "\n",
    "# Slave (pass 1): offset by the baseline\n",
    "pos1 = torch.zeros([nsweeps, 3], dtype=torch.float32, device=device)\n",
    "pos1[:, 0] = cross_track\n",
    "pos1[:, 1] = y_coords.to(device) + along_track\n",
    "pos1[:, 2] = altitude + dh\n",
    "\n",
    "# Aperture (azimuth) and range windows\n",
    "wa = torch.tensor(signal.get_window((\"taylor\", 3, 30), nsweeps, fftbins=False),\n",
    "                  dtype=torch.float32, device=device)\n",
    "wa /= wa.mean()\n",
    "wr = torch.tensor(signal.get_window((\"taylor\", 3, 30), nsamples, fftbins=False),\n",
    "                  dtype=torch.float32, device=device)\n",
    "wr /= wr.mean()\n",
    "\n",
    "# Window-weighted antenna phase centres\n",
    "apc0 = (wa[:, None] * pos0).mean(0) / wa.mean()\n",
    "apc1 = (wa[:, None] * pos1).mean(0) / wa.mean()\n",
    "baseline = (apc1 - apc0)\n",
    "print(\"baseline (m):\", baseline.cpu().numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Forward projection and range compression\n",
    "\n",
    "We simulate the slave's raw FMCW data with `projection_cart_2d_nufft`, then range-compress it (windowed, oversampled inverse FFT). The `data_fmod` term shifts the range spectrum so that, after the oversampled transform, it sits at DC, which minimises later interpolation error. This is the same image-formation recipe used in the interferometry example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "d0 = 0\n",
    "oversample = 3\n",
    "n = nsamples * oversample\n",
    "r_res = 3e8 / (2 * bw * oversample)\n",
    "# data_fmod centers the input data spectrum decreasing interpolation errors in backprojection\n",
    "data_fmod = -torch.pi * (1 - (oversample - 1) / oversample)\n",
    "\n",
    "# Forward project the slave pass to raw FMCW data (stop-and-go NUFFT)\n",
    "data = torchbp.ops.projection_cart_2d_nufft(\n",
    "    img, pos1, grid_proj, fc, fs, bw / tsweep, nsamples, d0, dem,\n",
    "    use_rvp=False, normalization=\"gamma\")[0]\n",
    "\n",
    "# Range compression: window + oversampled IFFT, then shift spectrum to DC\n",
    "w = wa[:, None] * wr[None, :]\n",
    "data = torch.fft.ifft(data * w, dim=-1, n=n)\n",
    "data = data * torch.exp(1j * data_fmod * torch.arange(data.shape[-1], device=device))[None, :]\n",
    "print(\"range-compressed data:\", tuple(data.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "## Forming the two reference images\n",
    "\n",
    "We backproject the same slave data twice onto the same polar grid, differing only in the origin:\n",
    "\n",
    "  * `slave_own`, centered on the slave's own phase centre (`apc1`). This is the image we will later interpolate.\n",
    "  * `slave_master`, centered on the master's phase centre (`apc0`). This is the ground truth that `polar_interp` must reproduce.\n",
    "\n",
    "Backprojection is done with `dealias=True`. A critically sampled polar image has a range spectrum that wraps around (aliases); de-aliasing unwraps it into a clean one-sided spectrum (see `torchbp.util.bp_polar_range_dealias`). After de-aliasing the spectrum occupies $[0,\\,2\\pi/\\texttt{signal\\_oversample}]$ rather than being centred at DC. Interpolation is most accurate on a DC-centred signal, so we pass `alias_fmod` to center the spectrum at DC. The same `alias_fmod` is later handed to `polar_interp` so it knows where the spectrum sits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "alias_fmod = float(-torch.pi / signal_oversample)\n",
    "print(f\"alias_fmod = {alias_fmod:.4f}\")\n",
    "\n",
    "def backproject(pos_world, center_xy):\n",
    "    # Center the grid under `center_xy` by shifting the platform horizontally.\n",
    "    origin = torch.tensor([[center_xy[0], center_xy[1], 0.0]], device=device)\n",
    "    pos_centered = pos_world - origin\n",
    "    return torchbp.ops.backprojection_polar_2d(\n",
    "        data, grid_polar, fc, r_res, pos_centered, d0=d0,\n",
    "        dealias=True, data_fmod=data_fmod, alias_fmod=alias_fmod)[0]\n",
    "\n",
    "t = time.time()\n",
    "slave_own = backproject(pos1, apc1[:2])     # centred on slave APC\n",
    "slave_master = backproject(pos1, apc0[:2])  # centred on master APC (ground truth)\n",
    "print(f\"Both images formed in {time.time() - t:.2f} s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "The two images are formed from identical data, so their magnitudes look the same. The difference lives entirely in the phase: each is referenced to a different origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "extent = [*grid_polar[\"r\"], *grid_polar[\"theta\"]]\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
    "for ax, sar, title in zip(axes, [slave_own, slave_master], [\"Centred on slave APC\", \"Centred on master APC\"]):\n",
    "    db = 20 * torch.log10(torch.abs(sar) + 1e-9).cpu().numpy()\n",
    "    m = db.max()\n",
    "    ax.imshow(db.T, origin=\"lower\", vmin=m - 30, vmax=m, extent=extent, aspect=\"auto\")\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(\"Range (m)\")\n",
    "axes[0].set_ylabel(\"Angle (sin radians)\")\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Image spectrum\n",
    "\n",
    "For best interpolation results, the image spectrum must not be aliased and it should be centered at zero. We used `dealias=True` and `alias_fmod` to guarantee this during image formation.\n",
    "\n",
    "Note how the spectrums of the images are slightly different due to baseline offset. `slave_own` positions are centered around image origin which leads to more compact spectrum than `slave_master` that has the baseline offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_own_db = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft2(slave_own)))).T\n",
    "m = np.max(spec_own_db)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
    "axes[0].imshow(spec_own_db, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 40, extent=[-0.5, 0.5, -0.5, 0.5])\n",
    "axes[0].set_xlabel(\"Range freq (cycles/sample)\")\n",
    "axes[0].set_ylabel(\"Azimuth freq (cycles/sample)\")\n",
    "axes[0].set_title(\"slave_own spectrum\");\n",
    "\n",
    "spec_master_db = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft2(slave_master)))).T\n",
    "\n",
    "axes[1].imshow(spec_master_db, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 40, extent=[-0.5, 0.5, -0.5, 0.5])\n",
    "axes[1].set_xlabel(\"Range freq (cycles/sample)\")\n",
    "axes[1].set_ylabel(\"Azimuth freq (cycles/sample)\")\n",
    "axes[1].set_title(\"slave_master spectrum\")\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "Without `dealias=True` and `alias_fmod` the spectrum of the image can be very badly aliased making it unsuitable to use for interpolation, although the amplitude of the image would be identical. `torchbp.util.bp_polar_range_alias` and `torchbp.util.bp_polar_range_dealias` can be used to alias and dealias already formed image.\n",
    "\n",
    "It's important for phase sensitive applications, such as interferometry, to restore the phase after interpolation using `torchbp.util.bp_polar_range_alias`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is equivalent to if BP was called with dealias=False and alias_fmod=0\n",
    "aliased = torchbp.util.bp_polar_range_alias(slave_own, apc1, fc=fc, grid_polar=grid_polar, alias_fmod=alias_fmod)\n",
    "\n",
    "spec_aliased_db = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft2(aliased)))).T\n",
    "m = np.max(spec_aliased_db)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
    "axes[0].imshow(spec_aliased_db, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 40, extent=[-0.5, 0.5, -0.5, 0.5])\n",
    "axes[0].set_xlabel(\"Range freq (cycles/sample)\")\n",
    "axes[0].set_ylabel(\"Azimuth freq (cycles/sample)\")\n",
    "axes[0].set_title(\"Aliased image spectrum\");\n",
    "\n",
    "db = 20 * torch.log10(torch.abs(aliased) + 1e-9).cpu().numpy()\n",
    "m = db.max()\n",
    "axes[1].imshow(db.T, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 30, extent=[*grid_polar[\"r\"], *grid_polar[\"theta\"]])\n",
    "axes[1].set_xlabel(\"Range (m)\")\n",
    "axes[1].set_ylabel(\"Azimuth (sin(degrees))\")\n",
    "axes[1].set_title(\"Aliased image amplitude\")\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "## Why re-centering is needed\n",
    "\n",
    "If we naively treat `slave_own` as if it were already on the master grid and take the phase difference against `slave_master`, we just get random noise. The two images disagree by the geometric phase of the baseline and the same pixel locations between to images don't corresponds to the same location in the scene. Overlaying them directly is meaningless."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "naive = torch.angle(slave_own * slave_master.conj()).cpu().numpy()\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.imshow(naive.T, origin=\"lower\", cmap=\"hsv\", extent=extent, aspect=\"auto\")\n",
    "plt.colorbar(label=\"Phase (rad)\")\n",
    "plt.xlabel(\"Range (m)\")\n",
    "plt.ylabel(\"Angle (sin radians)\")\n",
    "plt.title(\"Naive phase difference (different origins)\")\n",
    "\n",
    "phase_std = torch.angle(slave_own * slave_master.conj())\n",
    "print(f\"Phase std {torch.std(phase_std).item():.2f} rad\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "## Re-centering with `polar_interp`\n",
    "\n",
    "Now interpolate `slave_own` onto the master centre. We work in the master's centered coordinate frame:\n",
    "\n",
    "  * `origin_old` is where the slave image is *currently* referenced, the slave phase centre, which sits at `baseline` in the master frame. Its $z$ is the slave's own acquisition height.\n",
    "  * `origin_new` is the target centre, the master grid origin at horizontal $(0, 0)$, keeping the slave's acquisition height.\n",
    "\n",
    "`polar_interp` shifts every pixel from `origin_old` to `origin_new`, resampling the geometry and updating the carrier phase using `fc`. We use Lanczos resampling and pass the `alias_fmod` used during backprojection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "origin_old = torch.tensor([baseline[0].item(), baseline[1].item(), apc1[2].item()], device=device)\n",
    "origin_new = torch.tensor([0.0, 0.0, apc1[2].item()], device=device)\n",
    "\n",
    "t = time.time()\n",
    "slave_interp = torchbp.ops.polar_interp(\n",
    "    slave_own, origin_old=origin_old, origin_new=origin_new,\n",
    "    grid_polar=grid_polar, fc=fc, rotation=0,\n",
    "    grid_polar_new=grid_polar, method=(\"lanczos\", 8),\n",
    "    alias_fmod=alias_fmod)\n",
    "# polar_interp may return extra batch dim; squeeze to 2D\n",
    "slave_interp = slave_interp.squeeze()\n",
    "print(f\"polar_interp in {time.time() - t:.3f} s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "## Comparison against ground truth\n",
    "\n",
    "The interpolated image should now match the directly-backprojected `slave_master`. It won't be exactly similar, since there will be some interpolation error, but it should be close. We mask off the borders (where the interpolation kernel runs out of support) and look at three things:\n",
    "\n",
    "  * amplitude ratio, should be ~1\n",
    "  * residual phase, $\\angle(\\,s_\\text{interp}\\, s_\\text{master}^{*})$, should be ~0\n",
    "  * coherence, should be ~1\n",
    "\n",
    "The main source of error in this case is the aliasing during initial backprojection. Linear interpolation of the input data is used, which causes slight spreading of the spectrum which aliases during interpolation. However, this error is not very large and can be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = (torch.abs(slave_interp) > 0) & (torch.abs(slave_master) > 0)\n",
    "b = 20\n",
    "mask[:b, :] = False; mask[-b:, :] = False\n",
    "mask[:, :b] = False; mask[:, -b:] = False\n",
    "\n",
    "phase_res = torch.angle(slave_interp * slave_master.conj())\n",
    "amp_ratio = torch.abs(slave_interp) / (torch.abs(slave_master) + 1e-10)\n",
    "coh = torchbp.ops.coherence.coherence_2d(slave_interp, slave_master, (5, 5)).squeeze()\n",
    "\n",
    "print(f\"residual phase : mean {phase_res[mask].mean():+.4f}  std {phase_res[mask].std():.4f} rad\")\n",
    "print(f\"amplitude ratio: mean {amp_ratio[mask].mean():.4f}  std {amp_ratio[mask].std():.4f}\")\n",
    "print(f\"coherence      : mean {coh[mask].mean():.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "\n",
    "db_m = 20 * torch.log10(torch.abs(slave_master) + 1e-9)\n",
    "db_i = 20 * torch.log10(torch.abs(slave_interp) + 1e-9)\n",
    "im0 = axes[0].imshow((db_i - db_m).cpu().numpy().T, origin=\"lower\", cmap=\"RdBu_r\",\n",
    "                     vmin=-3, vmax=3, extent=extent, aspect=\"auto\")\n",
    "axes[0].set_title(\"Amplitude difference (dB)\")\n",
    "plt.colorbar(im0, ax=axes[0])\n",
    "\n",
    "im1 = axes[1].imshow(phase_res.cpu().numpy().T, origin=\"lower\", cmap=\"RdBu_r\",\n",
    "                     vmin=-0.5, vmax=0.5, extent=extent, aspect=\"auto\")\n",
    "axes[1].set_title(\"Residual phase (rad)\")\n",
    "plt.colorbar(im1, ax=axes[1])\n",
    "\n",
    "im2 = axes[2].imshow(coh.cpu().numpy().T, origin=\"lower\", vmin=0, vmax=1,\n",
    "                     extent=extent, aspect=\"auto\")\n",
    "axes[2].set_title(\"Coherence\")\n",
    "plt.colorbar(im2, ax=axes[2])\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlabel(\"Range (m)\")\n",
    "axes[0].set_ylabel(\"Angle (sin radians)\")\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29",
   "metadata": {},
   "source": [
    "## Spectrum of the interpolated image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_interp_db = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft2(slave_interp)))).T\n",
    "m = np.max(spec_interp_db)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
    "axes[0].imshow(spec_interp_db, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 40, extent=[-0.5, 0.5, -0.5, 0.5])\n",
    "axes[0].set_xlabel(\"Range freq (cycles/sample)\")\n",
    "axes[0].set_ylabel(\"Azimuth freq (cycles/sample)\")\n",
    "axes[0].set_title(\"Interpolated image spectrum\");\n",
    "\n",
    "axes[1].imshow(spec_master_db, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 40, extent=[-0.5, 0.5, -0.5, 0.5])\n",
    "axes[1].set_xlabel(\"Range freq (cycles/sample)\")\n",
    "axes[1].set_ylabel(\"Azimuth freq (cycles/sample)\")\n",
    "axes[1].set_title(\"Backprojected image spectrum\");\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31",
   "metadata": {},
   "source": [
    "## Interpolation kernel\n",
    "\n",
    "Linear interpolation is much faster, but it also has more interpolation error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = time.time()\n",
    "slave_interp = torchbp.ops.polar_interp(\n",
    "    slave_own, origin_old=origin_old, origin_new=origin_new,\n",
    "    grid_polar=grid_polar, fc=fc, rotation=0,\n",
    "    grid_polar_new=grid_polar, method=\"linear\",\n",
    "    alias_fmod=alias_fmod)\n",
    "# polar_interp may return a tuple / extra batch dim; squeeze to 2D\n",
    "if isinstance(slave_interp, tuple):\n",
    "    slave_interp = slave_interp[0]\n",
    "while slave_interp.dim() > 2:\n",
    "    slave_interp = slave_interp[0]\n",
    "print(f\"polar_interp in {time.time() - t:.3f} s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = (torch.abs(slave_interp) > 0) & (torch.abs(slave_master) > 0)\n",
    "b = 20\n",
    "mask[:b, :] = False; mask[-b:, :] = False\n",
    "mask[:, :b] = False; mask[:, -b:] = False\n",
    "\n",
    "phase_res = torch.angle(slave_interp * slave_master.conj())\n",
    "amp_ratio = torch.abs(slave_interp) / (torch.abs(slave_master) + 1e-10)\n",
    "coh = torchbp.ops.coherence.coherence_2d(slave_interp, slave_master, (5, 5)).squeeze()\n",
    "\n",
    "print(f\"residual phase : mean {phase_res[mask].mean():+.4f}  std {phase_res[mask].std():.4f} rad\")\n",
    "print(f\"amplitude ratio: mean {amp_ratio[mask].mean():.4f}  std {amp_ratio[mask].std():.4f}\")\n",
    "print(f\"coherence      : mean {coh[mask].mean():.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "\n",
    "db_m = 20 * torch.log10(torch.abs(slave_master) + 1e-9)\n",
    "db_i = 20 * torch.log10(torch.abs(slave_interp) + 1e-9)\n",
    "im0 = axes[0].imshow((db_i - db_m).cpu().numpy().T, origin=\"lower\", cmap=\"RdBu_r\",\n",
    "                     vmin=-3, vmax=3, extent=extent, aspect=\"auto\")\n",
    "axes[0].set_title(\"Amplitude difference (dB)\")\n",
    "plt.colorbar(im0, ax=axes[0])\n",
    "\n",
    "im1 = axes[1].imshow(phase_res.cpu().numpy().T, origin=\"lower\", cmap=\"RdBu_r\",\n",
    "                     vmin=-0.5, vmax=0.5, extent=extent, aspect=\"auto\")\n",
    "axes[1].set_title(\"Residual phase (rad)\")\n",
    "plt.colorbar(im1, ax=axes[1])\n",
    "\n",
    "im2 = axes[2].imshow(coh.cpu().numpy().T, origin=\"lower\", vmin=0, vmax=1,\n",
    "                     extent=extent, aspect=\"auto\")\n",
    "axes[2].set_title(\"Coherence\")\n",
    "plt.colorbar(im2, ax=axes[2])\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlabel(\"Range (m)\")\n",
    "axes[0].set_ylabel(\"Angle (sin radians)\")\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_interp_db = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft2(slave_interp)))).T\n",
    "m = np.max(spec_interp_db)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
    "axes[0].imshow(spec_interp_db, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 40, extent=[-0.5, 0.5, -0.5, 0.5])\n",
    "axes[0].set_xlabel(\"Range freq (cycles/sample)\")\n",
    "axes[0].set_ylabel(\"Azimuth freq (cycles/sample)\")\n",
    "axes[0].set_title(\"Linear interpolated image spectrum\");\n",
    "\n",
    "axes[1].imshow(spec_master_db, origin=\"lower\", aspect=\"auto\", vmax=m, vmin=m - 40, extent=[-0.5, 0.5, -0.5, 0.5])\n",
    "axes[1].set_xlabel(\"Range freq (cycles/sample)\")\n",
    "axes[1].set_ylabel(\"Azimuth freq (cycles/sample)\")\n",
    "axes[1].set_title(\"Backprojected image spectrum\");\n",
    "plt.tight_layout();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}