{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Interferometry with synthetic data\n",
    "\n",
    "Synthetic Aperture Radar interferometry (InSAR) recovers terrain height by comparing the phases of two SAR images acquired from slightly different positions. The phase difference between the two images (interferogram) depends on the path-length difference from each antenna to the ground, which depends on the local terrain height and the baseline separation of the two acquisitions.\n",
    "\n",
    "This example:\n",
    "  1. Builds a synthetic scene of random scatterers on top of a smooth digital elevation model (DEM).\n",
    "  2. Simulates two FMCW radar passes with a small vertical baseline using forward projection operator in `torchbp`.\n",
    "  3. Forms SAR images of both passes and computes their interferogram.\n",
    "  4. Compares the measured interferogram against an analytic prediction derived directly from the known geometry.\n",
    "  \n",
    "A complete InSAR processing chain would then unwrap the interferometricphase (which is only known modulo $2\\pi$) and convert it to height. `torchbp`does not include a phase-unwrapping algorithm. An external tool such as [SNAPHU](https://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu/) can be used for that step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "import torch.nn.functional as F\n",
    "import torchbp\n",
    "import torchbp.interferometry as ti\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import hamming\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",
    "We model a small FMCW radar flown along the $y$ axis at 100 m altitude. The two passes are separated by a 1 m vertical baseline (the second pass flies 1 m higher). The platform looks in the $+x$ direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 6e9            # RF center frequency (Hz)\n",
    "bw = 200e6          # Sweep bandwidth (Hz)\n",
    "tsweep = 100e-6     # Sweep length (s)\n",
    "fs = 20e6           # 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 = 100.0    # Platform altitude (m)\n",
    "baseline_z = 1.0    # Vertical baseline between the two passes (m)\n",
    "element_spacing = 0.25 * wl   # Azimuth sampling along the track"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "Polar imaging grid. `theta` is the sine of the azimuth angle. The number of range and azimuth bins is chosen from the range resolution and the synthetic aperture so the grid is slightly oversampled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "r0, r1 = 100.0, 250.0     # Ground-range extent (m)\n",
    "theta_limit = 0.4         # Azimuth half-extent (sine of angle)\n",
    "\n",
    "res = 3e8 / (2 * bw)      # Slant-range resolution\n",
    "res_interp = 1.3          # 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",
    "print(\"nr:\", nr, \"ntheta:\", ntheta, \"nsamples:\", nsamples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Scene\n",
    "\n",
    "Random scatterers over a smooth DEM. The scene is defined on a Cartesian grid on the ground ($z=0$) plane. Each grid cell is a scatterer with a random complex amplitude. Each scatterer's Z-coordinate is sampled from DEM so that they are at slightly different elevations. Because the DEM is smooth, the terrain height changes slowly across the scene and produces the gently curving interferometric fringes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.manual_seed(0)\n",
    "\n",
    "grid_proj = {\"x\": (60.0, 270.0), \"y\": (-100.0, 100.0), \"nx\": 384, \"ny\": 384}\n",
    "\n",
    "# Random complex amplitude per scatterer\n",
    "img = torch.randn([grid_proj[\"nx\"], grid_proj[\"ny\"]], dtype=torch.complex64, device=device)\n",
    "\n",
    "# Smooth, low-frequency DEM (height in metres)\n",
    "xg = torch.linspace(grid_proj[\"x\"][0], grid_proj[\"x\"][1], grid_proj[\"nx\"], device=device)\n",
    "yg = torch.linspace(grid_proj[\"y\"][0], grid_proj[\"y\"][1], grid_proj[\"ny\"], device=device)\n",
    "X, Y = torch.meshgrid(xg, yg, indexing=\"ij\")\n",
    "dem = (4.0 * torch.sin(2 * np.pi * (X - grid_proj[\"x\"][0]) / 180.0)\n",
    "       + 2.0 * torch.cos(2 * np.pi * Y / 140.0)).to(torch.float32)\n",
    "\n",
    "extent_proj = [*grid_proj[\"x\"], *grid_proj[\"y\"]]\n",
    "plt.figure()\n",
    "plt.imshow(dem.cpu().numpy().T, origin=\"lower\", extent=extent_proj, aspect=\"auto\")\n",
    "plt.colorbar(label=\"Height (m)\")\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"y (m)\")\n",
    "plt.title(\"Low-frequency DEM\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Platform trajectories\n",
    "\n",
    "Both passes fly straight along $y$ at $x=0$. The only difference is altitude, the second pass is `baseline_z` metres higher. Because both tracks share the same horizontal position, the two SAR images land on the same polar grid and no coregistration is needed. In practice the images would need to be coregistered since interferometry is very sensitive to position errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_coords = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps) * element_spacing\n",
    "pos0 = torch.zeros([nsweeps, 3], dtype=torch.float32, device=device)\n",
    "pos0[:, 1] = y_coords\n",
    "pos0[:, 2] = altitude\n",
    "\n",
    "pos1 = pos0.clone()\n",
    "pos1[:, 2] = altitude + baseline_z\n",
    "\n",
    "# Mean antenna phase-centre of each pass (the imaging origin)\n",
    "origin0 = pos0.mean(0)\n",
    "origin1 = pos1.mean(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Forward projection and image formation\n",
    "\n",
    "We simulate the raw radar data with `projection_cart_2d_nufft`, which generates radar measurement data from the scene using non-uniform FFT. There are two forward projectors in `torchbp`: `projection_cart_2d_nufft` is fast but uses the stop-and-go approximation, the platform is assumed to be stationary during each sweep. This is accurate for typical geometries and is what we use here. `projection_cart_2d` is a direct evaluation that also supports intra-sweep motion (`vel`). It is more general but much slower (typically about 10 to 100 times).\n",
    "\n",
    "After projection, the raw data is range-compressed (FFT with oversampling) and backprojected onto the polar grid. The `data_fmod` term centers the range spectrum around DC to minimize interpolation error. We backproject with `dealias=False` to not lose the interferometric phase (see `torchbp.util.bp_polar_range_dealias` for more info). Both passes are imaged on the same flat$z=0$ grid. Note that unlike slant range image formation, flat-plane backprojection interferogram does not include flat earth phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "oversample = 2\n",
    "data_fmod = -torch.pi * (1 - (oversample - 1) / oversample)\n",
    "n = nsamples * oversample\n",
    "r_res = 3e8 / (2 * bw * oversample)\n",
    "\n",
    "def make_image(pos):\n",
    "    # Forward project the scene to raw FMCW data (stop-and-go NUFFT)\n",
    "    data = torchbp.ops.projection_cart_2d_nufft(\n",
    "        img, pos, grid_proj, fc, fs, bw / tsweep, nsamples,\n",
    "        dem=dem, use_rvp=False, normalization=\"gamma\")[0]\n",
    "    # Range compression: window + oversampled IFFT\n",
    "    w = torch.tensor(hamming(data.shape[-1])[None, :], dtype=torch.float32, device=device)\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",
    "    # Backproject to the polar grid using FFBP, referenced to the flat z=0 plane\n",
    "    return torchbp.ops.backprojection_polar_2d(\n",
    "        data, grid_polar, fc, r_res, pos, dealias=False, data_fmod=data_fmod)[0]\n",
    "\n",
    "t = time.time()\n",
    "sar0 = make_image(pos0)\n",
    "sar1 = make_image(pos1)\n",
    "print(f\"Both passes formed in {time.time() - t:.2f} s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "The two SAR images look essentially identical in magnitude. The height information is hidden in their relative phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "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, [sar0, sar1], [\"Pass 0\", \"Pass 1\"]):\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": "16",
   "metadata": {},
   "source": [
    "Coherence of the images calcules how well the match to each other around a small patch. Coherence of 1 means they are identical, 0 means there is no correlation. Although simulation is ideal `torchbp.ops.coherence.coherence_2d` is not 1 everywhere because of the DEM causes phase difference in images. For proper coherence calculation DEM phase should be removed first.\n",
    "\n",
    "There's also `torchbp.ops.coherence.power_coherence_2d` that calculates coherence without phase and is not sensitive to DEM. However, because it doesn't use the phase it's lower quality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
    "\n",
    "im0 = axes[0].imshow(torchbp.ops.coherence.coherence_2d(sar0, sar1, (5,5)).cpu().numpy().T, vmin=0, vmax=1, origin=\"lower\", aspect=\"auto\", extent=extent)\n",
    "axes[0].set_xlabel(\"Range (m)\")\n",
    "axes[0].set_ylabel(\"Angle (sin radians)\")\n",
    "axes[0].set_title(\"Coherence\")\n",
    "plt.colorbar(im0, ax=axes[0])\n",
    "\n",
    "im1 = axes[1].imshow(torchbp.ops.coherence.power_coherence_2d(sar0, sar1, (5,5)).cpu().numpy().T, vmin=0, vmax=1, origin=\"lower\", aspect=\"auto\", extent=extent)\n",
    "axes[1].set_xlabel(\"Range (m)\")\n",
    "axes[1].set_ylabel(\"Angle (sin radians)\")\n",
    "axes[1].set_title(\"Power coherence\");\n",
    "plt.colorbar(im0, ax=axes[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Interferogram\n",
    "\n",
    "The interferogram is the product $s_0 \\, s_1^{*}$. Because the images were backprojected onto the flat $z=0$ plane, the flat-earth component is already removed. The interferometric phase is topographic and, to first order, proportional to terrain height. The fringes follow the shape of the DEM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "inf = sar0 * sar1.conj()\n",
    "inf_phase = torch.angle(inf).cpu().numpy()\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.imshow(inf_phase.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(\"Measured interferogram\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "### Phase unwrapping and DEM generation (not included)\n",
    "\n",
    "The interferometric phase above is wrapped, it is only known modulo $2\\pi$. Turning it into absolute height requires phase unwrapping, which is a separate, non-trivial step.`torchbp` does not provide a phase-unwrapping algorithm. In a full pipeline you would export the wrapped interferogram (and a coherence map) to a dedicated unwrapper such as SNAPHU, then convert the unwrapped phase to height. For the conversion `torchbp` has `torchbp.interferometry.phase_to_elevation_polar`, which divides the unwrapped phase by the local height sensitivity$\\tfrac{4\\pi}{\\lambda}(z_2/r_2 - z_1/r_1)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "## Analytic interferogram\n",
    "\n",
    "Because the geometry and DEM are known exactly, we can predict the interferogram directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_arr = torch.linspace(r0, r1, nr, device=device)\n",
    "theta_arr = torch.linspace(-theta_limit, theta_limit, ntheta, device=device)\n",
    "# pixel x, y ground coordinates\n",
    "x_w = r_arr[:, None] * torch.sqrt(1 - theta_arr[None, :] ** 2)\n",
    "y_w = r_arr[:, None] * theta_arr[None, :]\n",
    "r_ground = r_arr[:, None].expand_as(x_w)\n",
    "\n",
    "def dem_at(xq, yq):\n",
    "    xn = 2 * (xq - grid_proj[\"x\"][0]) / (grid_proj[\"x\"][1] - grid_proj[\"x\"][0]) - 1\n",
    "    yn = 2 * (yq - grid_proj[\"y\"][0]) / (grid_proj[\"y\"][1] - grid_proj[\"y\"][0]) - 1\n",
    "    gc = torch.stack([yn, xn], dim=-1)[None]\n",
    "    return F.grid_sample(dem[None, None], gc, mode=\"bilinear\", align_corners=True)[0, 0]\n",
    "\n",
    "R0 = torch.sqrt(r_ground ** 2 + altitude ** 2)\n",
    "# Solve for the scatterer's range in slave image.\n",
    "# Layover/foreshortening shifts it in ground range.\n",
    "rho = r_ground.clone()\n",
    "for _ in range(6):\n",
    "    s = rho / r_ground # radial scale toward the true position\n",
    "    hs = dem_at(s * x_w, s * y_w)\n",
    "    rho = torch.sqrt(torch.clamp(R0 ** 2 - (altitude - hs) ** 2, min=0))\n",
    "R1 = torch.sqrt(rho ** 2 + (altitude + baseline_z - hs) ** 2)\n",
    "\n",
    "phase_full = (4 * np.pi / wl) * (R1 - R0) # total geometric phase\n",
    "phase_flat = (4 * np.pi / wl) * (torch.sqrt(r_ground ** 2 + (altitude + baseline_z) ** 2) - R0) # flat earth\n",
    "analytic = phase_full - phase_flat # topographic phase\n",
    "analytic_wrapped = torch.angle(torch.exp(1j * analytic))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "## Comparison\n",
    "\n",
    "The measured and analytic interferograms show the same topographic fringe pattern. After flattening with the analytic phase, the residual has still slight error due to layover. For more accurate DEM generation the images should be formed on a DEM. If no DEM exists, then first pass can be made without DEM and then the solved DEM is used for second pass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = inf_phase - analytic.cpu().numpy()\n",
    "residual = np.angle(np.exp(1j * diff))\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "axes[0].imshow(inf_phase.T, origin=\"lower\", cmap=\"hsv\", extent=extent, aspect=\"auto\")\n",
    "axes[0].set_title(\"Measured\")\n",
    "axes[1].imshow(analytic_wrapped.cpu().numpy().T, origin=\"lower\", cmap=\"hsv\", extent=extent, aspect=\"auto\")\n",
    "axes[1].set_title(\"Analytic\")\n",
    "im = axes[2].imshow(residual.T, origin=\"lower\", cmap=\"hsv\", vmin=-np.pi, vmax=np.pi,\n",
    "                    extent=extent, aspect=\"auto\")\n",
    "axes[2].set_title(\"Residual\")\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": "25",
   "metadata": {},
   "source": [
    "Speckle makes the interferogram noisy. It can filtered with multilooking (spatial averaging) or Goldstein phase filter. Both trade resolution for less speckle noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "inf_gold = torchbp.interferometry.goldstein_filter(inf)\n",
    "inf_multilook = torchbp.ops.multilook_polar(inf, (2,2), grid_polar)[0]\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "axes[0].imshow((torch.angle(inf_gold)).T.cpu().numpy(), origin=\"lower\", cmap=\"hsv\", extent=extent, aspect=\"auto\")\n",
    "axes[0].set_xlabel(\"Range (m)\")\n",
    "axes[0].set_ylabel(\"Angle (sin radians)\")\n",
    "axes[0].set_title(\"Goldstein filtered\")\n",
    "\n",
    "axes[1].imshow((torch.angle(inf_multilook)).T.cpu().numpy(), origin=\"lower\", cmap=\"hsv\", extent=extent, aspect=\"auto\")\n",
    "axes[1].set_xlabel(\"Range (m)\")\n",
    "axes[1].set_ylabel(\"Angle (sin radians)\")\n",
    "axes[1].set_title(\"Multilooked\")\n",
    "\n",
    "axes[2].imshow(analytic_wrapped.cpu().numpy().T, origin=\"lower\", cmap=\"hsv\", extent=extent, aspect=\"auto\")\n",
    "axes[2].set_title(\"Analytic\")\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27",
   "metadata": {},
   "source": [
    "## Cartesian projection\n",
    "\n",
    "Above images are in pseudo-polar coordinates. We can project them to Cartesian plane to compare with DEM. Interferograms are still wrapped, but their shape can be compared to DEM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "analytic_cart = torchbp.ops.polar_to_cart(analytic_wrapped, origin0, grid_polar, grid_proj, fc=0)[0]\n",
    "gold_cart = torchbp.ops.polar_to_cart(torch.angle(inf_gold), origin0, grid_polar, grid_proj, fc=0)[0]\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "axes[0].imshow(analytic_cart.T.cpu().numpy(), origin=\"lower\", extent=extent_proj, aspect=\"auto\", vmin=-np.pi, vmax=np.pi)\n",
    "axes[0].set_xlabel(\"X (m)\")\n",
    "axes[0].set_ylabel(\"Y (m)\")\n",
    "axes[0].set_title(\"Analytic interferogram\")\n",
    "\n",
    "axes[1].imshow(gold_cart.T.cpu().numpy(), origin=\"lower\", extent=extent_proj, aspect=\"auto\", vmin=-np.pi, vmax=np.pi)\n",
    "axes[1].set_xlabel(\"X (m)\")\n",
    "axes[1].set_ylabel(\"Y (m)\")\n",
    "axes[1].set_title(\"Goldstein filtered interferogram\")\n",
    "\n",
    "axes[2].imshow(dem.cpu().numpy().T, origin=\"lower\", extent=extent_proj, aspect=\"auto\")\n",
    "axes[2].set_xlabel(\"X (m)\")\n",
    "axes[2].set_ylabel(\"Y (m)\")\n",
    "axes[2].set_title(\"DEM\");"
   ]
  }
 ],
 "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
}