{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Cartesian factorized backprojection (CFBP)\n",
    "\n",
    "Like FFBP (see the ffbp example), CFBP reduces the $\\mathcal{O}(N^2)$ cost of direct\n",
    "backprojection by splitting the aperture into subapertures, imaging each one on a\n",
    "coarse grid, and coherently merging the images in a binary tree. The difference is the\n",
    "grid the factorization runs on:\n",
    "\n",
    "- **FFBP** forms each subaperture image on its own local polar grid, centered on the\n",
    "  subaperture. Merging requires interpolating between shifted polar frames, a 2D\n",
    "  interpolation with a finite kernel (`ffbp_merge2`).\n",
    "- **CFBP** forms every subaperture image on the same global Cartesian grid as the\n",
    "  final image, only coarsely sampled in the cross-range direction y. A subaperture\n",
    "  image is not bandlimited in y as-is because of the carrier phase\n",
    "  $\\exp(j 4\\pi f_c d / c)$, but the carrier is a known function of geometry, so it can\n",
    "  be removed by multiplying with a reference computed from the subaperture center.\n",
    "  The demodulated image is bandlimited, so the coarse y sampling is enough, and since\n",
    "  all subimages share the same grid extent, merging only needs interpolation along y-axis\n",
    "  plus a phase re-reference.\n",
    "\n",
    "The library exposes this as a single call, `torchbp.ops.cfbp`. In this notebook the\n",
    "algorithm is built by hand from `backprojection_cart_2d` and a few lines of torch to\n",
    "demonstrate how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import hamming\n",
    "\n",
    "import torchbp\n",
    "from torchbp.util import center_pos\n",
    "from torchbp.ops import backprojection_cart_2d\n",
    "\n",
    "device = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n",
    "print(\"Device:\", device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Scene and radar parameters\n",
    "\n",
    "A single point target imaged on a heavily oversampled Cartesian grid zoomed around the\n",
    "target so the point-spread function (PSF) is smooth and easy to inspect. Same radar\n",
    "parameters as the FFBP example. The platform track is along y (cross-range), boresight\n",
    "along x (range)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 6e9          # RF center frequency (Hz)\n",
    "bw = 200e6        # RF bandwidth (Hz), positive for rising sweep\n",
    "tsweep = 100e-6   # Sweep length (s)\n",
    "fs = 2e6          # Sampling frequency (Hz)\n",
    "nsweeps = 512     # Number of pulses in the aperture\n",
    "oversample = 2    # Range FFT oversampling (reduces interpolation error)\n",
    "nsamples = int(fs * tsweep)\n",
    "\n",
    "# Heavily oversampled Cartesian grid, zoomed on the target.\n",
    "nx, ny = 256, 512\n",
    "grid = {\"x\": (96.0, 104.0), \"y\": (-4.0, 4.0), \"nx\": nx, \"ny\": ny}\n",
    "\n",
    "# Single point target at 100 m, broadside.\n",
    "target_pos = torch.tensor([[100.0, 0.0, 0.0]], dtype=torch.float32, device=device)\n",
    "target_rcs = torch.ones((1, 1), dtype=torch.float32, device=device)\n",
    "\n",
    "# Straight platform track along y at 50 m height. Element spacing = lambda/4.\n",
    "pos = torch.zeros((nsweeps, 3), dtype=torch.float32, device=device)\n",
    "pos[:, 1] = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps, device=device) * 0.25 * 3e8 / fc\n",
    "pos[:, 2] = 50.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Simulate and range-compress the data\n",
    "\n",
    "Identical to the backprojection example. Generate raw data, apply a range window, and\n",
    "FFT to range-compress. `data_fmod` shifts the data spectrum to DC to reduce the linear\n",
    "interpolation error inside backprojection. Unlike FFBP there is no `alias_fmod`.\n",
    "The subaperture images are fully demodulated in 2D by the carrier reference, so no\n",
    "separate range-spectrum bookkeeping is needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw = torchbp.util.generate_fmcw_data(target_pos, target_rcs, pos, fc, bw, tsweep, fs)\n",
    "w = torch.tensor(hamming(raw.shape[-1])[None, :], dtype=torch.float32, device=device)\n",
    "data = torch.fft.ifft(raw * w, dim=-1, n=nsamples * oversample)\n",
    "\n",
    "data_fmod = -math.pi * (1 - (oversample - 1) / oversample)\n",
    "data = data * torch.exp(1j * data_fmod * torch.arange(data.shape[-1], device=device))[None, :]\n",
    "\n",
    "r_res = 3e8 / (2 * abs(bw) * oversample)  # range bin size in the data\n",
    "print(f\"r_res={r_res:.3f} m, data_fmod={data_fmod:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "Helpers: display a Cartesian image in dB, and the carrier reference. The image pixel\n",
    "phase is $\\pi k_\\text{eff} d$ with $k_\\text{eff} = 4 f_c / c - f_\\text{mod} / (\\pi\n",
    "\\Delta_r)$ (the `f_mod` term undoes phase of `data_fmod` term) and $d$ the distance from the antenna phase center (APC) to the\n",
    "pixel. `carrier_ref` evaluates $\\exp(\\pm j \\pi k_\\text{eff} d)$ on a grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "C0 = 299792458.0\n",
    "keff = 4 * fc / C0 - data_fmod / (math.pi * r_res)\n",
    "\n",
    "def grid_distance(g, origin, z, dtype=torch.float32):\n",
    "    x0, x1 = g[\"x\"]; y0, y1 = g[\"y\"]\n",
    "    dx = (x1 - x0) / g[\"nx\"]; dy = (y1 - y0) / g[\"ny\"]\n",
    "    # 64-bit precision for more accurate phase computation\n",
    "    x = x0 - float(origin[0]) + dx * torch.arange(g[\"nx\"], device=device, dtype=torch.float64)\n",
    "    y = y0 - float(origin[1]) + dy * torch.arange(g[\"ny\"], device=device, dtype=torch.float64)\n",
    "    d = torch.sqrt(x[:, None] ** 2 + y[None, :] ** 2 + float(z) ** 2)\n",
    "    return d.to(dtype)\n",
    "\n",
    "def carrier_ref(g, origin, z, sign):\n",
    "    d = grid_distance(g, origin, z, dtype=torch.float64)\n",
    "    ph = torch.remainder(sign * keff * d, 2.0).float() * torch.pi\n",
    "    return torch.polar(torch.ones_like(ph), ph)\n",
    "\n",
    "extent = [*grid[\"y\"], *grid[\"x\"]]\n",
    "\n",
    "def show_cart(img, title, ax=None, dyn=30, ref_max=None):\n",
    "    img = img.detach().cpu()\n",
    "    db = 20 * torch.log10(torch.abs(img) + 1e-12)\n",
    "    m = db.max() if ref_max is None else ref_max\n",
    "    if ax is None:\n",
    "        _, ax = plt.subplots()\n",
    "    im = ax.imshow(db.numpy(), origin=\"lower\", vmin=m - dyn, vmax=m,\n",
    "                   extent=extent, aspect=\"auto\")\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(\"Cross-range y (m)\")\n",
    "    ax.set_ylabel(\"Range x (m)\")\n",
    "    return m\n",
    "\n",
    "def relerr(a, b):\n",
    "    return (torch.linalg.norm(a - b) / torch.linalg.norm(a)).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Reference: Direct backprojection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_bp = backprojection_cart_2d(data, grid, fc, r_res, pos, data_fmod=data_fmod).squeeze(0)\n",
    "\n",
    "m = show_cart(img_bp, \"Direct backprojection (reference)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Why demodulation is needed\n",
    "\n",
    "The image carrier has the local y frequency $k_\\text{eff} (y - y_\\text{APC}) /\n",
    "(2 d)$ cycles/m at each target, so on a wide swath the carrier spreads the image\n",
    "spectrum far beyond what the coarse subaperture grid can sample. The zoomed\n",
    "single-target grid of this notebook is too small (and too oversampled) to show this,\n",
    "so this cell uses a separate wide demo scene with point targets across the swath.\n",
    "A half-aperture image is formed on a fine grid and its spectrum along y is plotted\n",
    "before and after removing the carrier reference: each target contributes a spectral\n",
    "line at its carrier frequency, and demodulation collapses all of them to baseband,\n",
    "leaving only the subaperture azimuth bandwidth. The dashed lines mark the Nyquist\n",
    "limit of a twice-coarser grid: without demodulation the wide-swath image would alias\n",
    "there, demodulated it fits with a large margin. This is what makes the coarse y\n",
    "sampling of the factorization valid, independent of the scene size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "demo_targets = torch.zeros((9, 3), dtype=torch.float32, device=device)\n",
    "demo_targets[:, 0] = 100.0\n",
    "demo_targets[:, 1] = torch.linspace(-36.0, 36.0, 9, device=device)\n",
    "demo_rcs = torch.ones((9, 1), dtype=torch.float32, device=device)\n",
    "raw_d = torchbp.util.generate_fmcw_data(demo_targets, demo_rcs, pos, fc, bw, tsweep, fs)\n",
    "data_d = torch.fft.ifft(raw_d * w, dim=-1, n=nsamples * oversample)\n",
    "data_d = data_d * torch.exp(1j * data_fmod * torch.arange(data_d.shape[-1], device=device))[None, :]\n",
    "\n",
    "ny_demo = 3072\n",
    "grid_demo = {\"x\": (96.0, 104.0), \"y\": (-40.0, 40.0), \"nx\": 128, \"ny\": ny_demo}\n",
    "\n",
    "p_half = pos[: nsweeps // 2]\n",
    "p_local, origin_row = center_pos(p_half)\n",
    "origin_half = origin_row[0]\n",
    "z_half = p_half[:, 2].mean()\n",
    "g_shift = dict(grid_demo,\n",
    "               x=(grid_demo[\"x\"][0] - float(origin_half[0]), grid_demo[\"x\"][1] - float(origin_half[0])),\n",
    "               y=(grid_demo[\"y\"][0] - float(origin_half[1]), grid_demo[\"y\"][1] - float(origin_half[1])))\n",
    "img_half = backprojection_cart_2d(data_d[: nsweeps // 2], g_shift, fc, r_res, p_local,\n",
    "                                  data_fmod=data_fmod).squeeze(0)\n",
    "img_half_demod = img_half * carrier_ref(grid_demo, origin_half, z_half, -1.0)\n",
    "\n",
    "dy_demo = (grid_demo[\"y\"][1] - grid_demo[\"y\"][0]) / ny_demo\n",
    "fy = np.fft.fftshift(np.fft.fftfreq(ny_demo, d=dy_demo))\n",
    "spec = lambda im: 20 * np.log10(np.abs(np.fft.fftshift(np.fft.fft(im.cpu().numpy(), axis=-1), axes=-1)).mean(axis=0) + 1e-9)\n",
    "\n",
    "s0 = spec(img_half); s1 = spec(img_half_demod)\n",
    "plt.plot(fy, s0 - s0.max(), label=\"with carrier\")\n",
    "plt.plot(fy, s1 - s0.max(), label=\"demodulated\")\n",
    "nyq_sub = 1 / (2 * dy_demo * 2)  # Nyquist of a 2x coarser grid\n",
    "plt.axvline(-nyq_sub, color=\"k\", ls=\"--\", lw=1); plt.axvline(nyq_sub, color=\"k\", ls=\"--\", lw=1)\n",
    "plt.xlim(-1.6 * nyq_sub, 1.6 * nyq_sub); plt.ylim(-60, 3)\n",
    "plt.xlabel(\"y frequency (cycles/m)\"); plt.ylabel(\"dB\")\n",
    "plt.title(\"Wide-swath subaperture image spectrum along y\"); plt.legend(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Step 1: Backproject the subapertures (coarse y)\n",
    "\n",
    "Split the aperture in two halves. Each half is backprojected with\n",
    "`backprojection_cart_2d` on the global grid extent with half the y samples, using\n",
    "positions centered on that half's mean APC (`center_pos`). The\n",
    "grid is shifted by the same origin, so pixels stay at the same global coordinates.\n",
    "Then the carrier referenced to the subaperture APC is removed. A small internal\n",
    "oversampling factor in y limits the interpolation error in the later merge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "OSR_Y = 1.4  # internal y oversampling of subaperture grids\n",
    "\n",
    "def backproject_subaperture(sl, ny_sub):\n",
    "    p = pos[sl]\n",
    "    p_local, origin_row = center_pos(p)\n",
    "    origin = origin_row[0]\n",
    "    z = p[:, 2].mean()\n",
    "    g = dict(grid, ny=ny_sub)\n",
    "    g_shift = dict(g,\n",
    "                   x=(g[\"x\"][0] - float(origin[0]), g[\"x\"][1] - float(origin[0])),\n",
    "                   y=(g[\"y\"][0] - float(origin[1]), g[\"y\"][1] - float(origin[1])))\n",
    "    img = backprojection_cart_2d(data[sl], g_shift, fc, r_res, p_local,\n",
    "                                 data_fmod=data_fmod).squeeze(0)\n",
    "    img = img * carrier_ref(g, origin, z, -1.0)  # demodulate\n",
    "    return {\"img\": img, \"grid\": g, \"origin\": origin, \"z\": z}\n",
    "\n",
    "half = nsweeps // 2\n",
    "subA = backproject_subaperture(slice(0, half), int(OSR_Y * ny // 2))\n",
    "subB = backproject_subaperture(slice(half, nsweeps), int(OSR_Y * ny // 2))\n",
    "\n",
    "print(\"subaperture image shape:\", tuple(subA[\"img\"].shape), \"(full grid is\", (nx, ny), \")\")\n",
    "fig, axs = plt.subplots(1, 2, figsize=(11, 4))\n",
    "show_cart(subA[\"img\"], \"Subaperture A (first half)\", ax=axs[0])\n",
    "show_cart(subB[\"img\"], \"Subaperture B (second half)\", ax=axs[1])\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "Each subaperture image is focused but the cross-range resolution is half of the\n",
    "reference since only half of the data is used. Unlike FFBP, both images are already in\n",
    "the same global coordinate frame, the target sits at the same (x, y) in both and only\n",
    "the demodulation references differ."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "## Step 2: Merge\n",
    "\n",
    "1. Upsample along y to the output grid with FFT zero-padding. Because all grids share\n",
    "   the same extent and samples sit at DFT sample positions, this interpolation is\n",
    "   exact for the bandlimited demodulated image. `torchbp.ops.cfbp` also supports faster\n",
    "   finite-length interpolation kernel.\n",
    "2. Re-reference the demodulation carrier from the subaperture APC to the merged APC by\n",
    "   multiplying with $\\exp(j \\pi k_\\text{eff} (d_\\text{sub} - d_\\text{merged}))$.\n",
    "3. Sum.\n",
    "\n",
    "After the final remodulation the result matches direct backprojection up to interpolation error.\n",
    "The upsampler below is the upsampling branch of the library's\n",
    "`torchbp.ops.cfbp._fft_resample_dim`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fft_upsample_y(img, ny_new):\n",
    "    ny_old = img.shape[-1]\n",
    "    if ny_new == ny_old:\n",
    "        return img\n",
    "    X = torch.fft.fft(img, dim=-1)\n",
    "    h = ny_old // 2\n",
    "    Y = X.new_zeros((*X.shape[:-1], ny_new))\n",
    "    if ny_old % 2 == 0:\n",
    "        Y[..., :h] = X[..., :h]\n",
    "        Y[..., h] = 0.5 * X[..., h]          # split the Nyquist bin\n",
    "        Y[..., ny_new - h] = 0.5 * X[..., h]\n",
    "        Y[..., ny_new - h + 1:] = X[..., h + 1:]\n",
    "    else:\n",
    "        Y[..., :h + 1] = X[..., :h + 1]\n",
    "        Y[..., ny_new - h:] = X[..., h + 1:]\n",
    "    return torch.fft.ifft(Y, dim=-1) * (ny_new / ny_old)\n",
    "\n",
    "def merge(a, b, ny_new):\n",
    "    new_origin = 0.5 * (a[\"origin\"] + b[\"origin\"])\n",
    "    new_z = 0.5 * (a[\"z\"] + b[\"z\"])\n",
    "    g_new = dict(grid, ny=ny_new)\n",
    "    d_new = grid_distance(g_new, new_origin, new_z)\n",
    "    out = 0\n",
    "    for s in (a, b):\n",
    "        img = fft_upsample_y(s[\"img\"], ny_new)\n",
    "        ph = (torch.pi * keff) * (grid_distance(g_new, s[\"origin\"], s[\"z\"]) - d_new)\n",
    "        out = out + img * torch.polar(torch.ones_like(ph), ph)\n",
    "    return {\"img\": out, \"grid\": g_new, \"origin\": new_origin, \"z\": new_z}\n",
    "\n",
    "merged = merge(subA, subB, ny)\n",
    "img_cfbp2 = merged[\"img\"] * carrier_ref(grid, merged[\"origin\"], merged[\"z\"], 1.0)  # remodulate\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(11, 4))\n",
    "show_cart(img_bp, \"Direct backprojection\", ax=axs[0], ref_max=m)\n",
    "show_cart(img_cfbp2, \"CFBP, 2 subapertures (manual)\", ax=axs[1], ref_max=m)\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "print(f\"relative error vs direct BP: {relerr(img_bp, img_cfbp2):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "The merged image matches the reference. The cross-range cut through the peak overlays the direct result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "pk = np.unravel_index(torch.abs(img_bp).cpu().argmax(), img_bp.shape)\n",
    "y_axis = np.linspace(grid[\"y\"][0], grid[\"y\"][1], ny)\n",
    "cut_bp = 20 * torch.log10(torch.abs(img_bp[pk[0]]).cpu() + 1e-12)\n",
    "cut_cf = 20 * torch.log10(torch.abs(img_cfbp2[pk[0]]).cpu() + 1e-12)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(12, 4))\n",
    "axs[0].plot(y_axis, cut_bp - cut_bp.max(), label=\"direct BP\")\n",
    "axs[0].plot(y_axis, cut_cf - cut_bp.max(), \"--\", label=\"CFBP (manual)\")\n",
    "axs[0].set_xlim(-1.5, 1.5); axs[0].set_ylim(-40, 1)\n",
    "axs[0].set_xlabel(\"Cross-range y (m)\"); axs[0].set_ylabel(\"dB\")\n",
    "axs[0].set_title(\"Cross-range cut through peak\"); axs[0].legend()\n",
    "show_cart(img_bp - img_cfbp2, \"Difference (CFBP - direct)\", ax=axs[1], ref_max=m, dyn=80)\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "## Recursion\n",
    "\n",
    "As with FFBP, the real speedup comes from recursing: four subapertures are\n",
    "backprojected on grids with a quarter of the y samples and merged in a two-level\n",
    "binary tree. Each level of the tree doubles the cross-range resolution and the number\n",
    "of y samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "NSUB = 4\n",
    "n = nsweeps // NSUB\n",
    "leaves = [backproject_subaperture(slice(i * n, (i + 1) * n), int(OSR_Y * ny // NSUB))\n",
    "          for i in range(NSUB)]\n",
    "\n",
    "# Level 1: merge pairs onto y-oversampled intermediate grids.\n",
    "mid0 = merge(leaves[0], leaves[1], int(OSR_Y * ny // 2))\n",
    "mid1 = merge(leaves[2], leaves[3], int(OSR_Y * ny // 2))\n",
    "\n",
    "# Level 2: merge the intermediates onto the full output grid and remodulate.\n",
    "merged4 = merge(mid0, mid1, ny)\n",
    "img_cfbp4 = merged4[\"img\"] * carrier_ref(grid, merged4[\"origin\"], merged4[\"z\"], 1.0)\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 4))\n",
    "show_cart(leaves[0][\"img\"], \"Leaf 0 (1/4 aperture)\", ax=axs[0])\n",
    "show_cart(mid0[\"img\"], \"After level-1 merge (1/2)\", ax=axs[1])\n",
    "show_cart(img_cfbp4, \"After level-2 merge (full)\", ax=axs[2], ref_max=m)\n",
    "plt.tight_layout(); plt.show()\n",
    "\n",
    "print(f\"4-subaperture CFBP, relative error vs direct BP: {relerr(img_bp, img_cfbp4):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "## Library op `torchbp.ops.cfbp`\n",
    "\n",
    "Everything above is what `torchbp.ops.cfbp` does internally, plus the y-axis guard band to reduce interpolation error at the edges and\n",
    "support for uneven splits and arbitrary grid sizes. `stages` is the number of recursion levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_lib1 = torchbp.ops.cfbp(data, grid, fc, r_res, pos, stages=1,\n",
    "                            data_fmod=data_fmod).squeeze(0)\n",
    "img_lib2 = torchbp.ops.cfbp(data, grid, fc, r_res, pos, stages=2,\n",
    "                            data_fmod=data_fmod).squeeze(0)\n",
    "\n",
    "print(f\"stages=1 vs direct BP: {relerr(img_bp, img_lib1):.4f}\")\n",
    "print(f\"stages=2 vs direct BP: {relerr(img_bp, img_lib2):.4f}\")\n",
    "print(f\"stages=1 vs manual 2-subaperture: {relerr(img_cfbp2, img_lib1):.4f}\")\n",
    "print(f\"stages=2 vs manual 4-subaperture: {relerr(img_cfbp4, img_lib2):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "## Range-adaptive variant `torchbp.ops.cfbp_adaptive`\n",
    "\n",
    "A single global `stages` is wasteful when the scene spans a wide range of ground\n",
    "ranges. The y-bandwidth of a demodulated subaperture grows with the angular extent of\n",
    "the aperture seen from a pixel, which is largest at near range, so the coarsest grid a\n",
    "deep recursion can use is set by the closest range while far ranges could tolerate much\n",
    "deeper factorization. Because the Cartesian scheme never resamples the x (range) axis,\n",
    "the scene can be split into ground-range blocks with no seam error, and each block run\n",
    "at its own y-density and stage count: `cfbp_adaptive` picks an integer y-oversampling\n",
    "multiplier `k(x)` and recursion depth `s(x)` per block, runs the `cfbp` tree on a grid\n",
    "with `k * ny` y-samples, and decimates by taking every k-th sample, which lands exactly\n",
    "on the requested output positions. Far ranges keep the full factorization speedup while\n",
    "near ranges automatically fall back to shallower recursion, so the whole image matches\n",
    "direct backprojection to normal CFBP accuracy. `stages` in this case is the maximum depth."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "## Benefit and downside\n",
    "\n",
    "Compared to direct backprojection, CFBP has the same $\\mathcal{O}(N \\log N)$ scaling\n",
    "as FFBP. Compared to FFBP it has two attractive properties: the merge interpolation is\n",
    "exact (FFT zero-padding on nested grids) rather than a finite-kernel approximation,\n",
    "and the output is directly on the Cartesian grid with no final polar-to-Cartesian\n",
    "resampling step.\n",
    "\n",
    "The downside is a geometric one. The demodulated subaperture images contain, in\n",
    "addition to the azimuth bandwidth that shrinks with the subaperture length, the range\n",
    "envelope bandwidth projected onto y, approximately $\\sin(\\theta_\\text{max}) B$ where\n",
    "$\\theta_\\text{max}$ is the largest scene angle from boresight and $B$ the range\n",
    "bandwidth of the data. This term does not shrink when the aperture is split, so wide\n",
    "angular extents combined with little range oversampling need a larger `oversample_y`\n",
    "or fewer `stages`. Polar FFBP does not have this limitation (the envelope varies along\n",
    "its range axis, which is never split), which makes it the better choice for very\n",
    "wide-angle imaging."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}