stitching_satellite_patches.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import cv2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage.registration import phase_cross_correlation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Before running the code below, put all satellite patches together and replace the directory below "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sat_image_root = './Satellite/' # path to the satellite patches\n",
    "_, _, sat_filenames = next(os.walk(sat_image_root))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sat_map = np.zeros([20000, 17000, 3])\n",
    "image1 = cv2.imread(sat_image_root + '619400_5736195.png')\n",
    "sat_map[600:1200, 600:1200, :] = image1\n",
    "y = 900\n",
    "x = 900\n",
    "for j in range(280):\n",
    "    print(j)\n",
    "    if j %2 == 0:\n",
    "        for i in range(0, 320):\n",
    "            if j == 0 and i == 0:\n",
    "                pass\n",
    "            else:\n",
    "                image2 = cv2.imread(sat_image_root + str(619400+int(j*5))+'_'+str(int(5736195-i*5))+'.png')\n",
    "                if image2 is not None:\n",
    "                    shift, _, _ = phase_cross_correlation(image1[:,:,0], image2[:,:,0])\n",
    "                    y_offset, x_offset = shift\n",
    "                    y += int(y_offset)\n",
    "                    x += int(x_offset)\n",
    "                    sat_map[y-300:y+300, x-300:x+300, :] = image2\n",
    "                    image1 = image2\n",
    "    else:\n",
    "        for i in range(319,-1,-1):\n",
    "            image2 = cv2.imread(sat_image_root + str(619400+int(j*5))+'_'+str(int(5736195-i*5))+'.png')\n",
    "            if image2 is not None:\n",
    "                shift, _, _ = phase_cross_correlation(image1[:,:,0], image2[:,:,0])\n",
    "                y_offset, x_offset = shift\n",
    "                y += int(y_offset)\n",
    "                x += int(x_offset)\n",
    "                sat_map[y-300:y+300, x-300:x+300, :] = image2\n",
    "                image1 = image2\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write the stitched satellite map to disk\n",
    "cv2.imwrite('./satellite_map_new.png', sat_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_bl = cv2.imread(sat_image_root+'619400_5734600.png') # bottom left\n",
    "image_br = cv2.imread(sat_image_root+'620795_5734600.png') # bottom right\n",
    "image_tr = cv2.imread(sat_image_root+'620795_5736195.png') # top right\n",
    "image_c = cv2.imread(sat_image_root+'620100_5735400.png') # center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_bl = sat_map[18168-300:18168+300, 492-300:492+300, :]\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(cv2.cvtColor(image_bl, cv2.COLOR_BGR2RGB))\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(cv2.cvtColor(map_bl.astype(np.uint8), cv2.COLOR_BGR2RGB))\n",
    "\n",
    "shift, _, _ = phase_cross_correlation(map_bl[:,:,0], image_bl[:,:,0])\n",
    "y_offset, x_offset = shift\n",
    "print(y_offset, x_offset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_br = sat_map[18528-300:18528+300, 15553-300:15553+300, :]\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(cv2.cvtColor(image_br, cv2.COLOR_BGR2RGB))\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(cv2.cvtColor(map_br.astype(np.uint8), cv2.COLOR_BGR2RGB))\n",
    "\n",
    "shift, _, _ = phase_cross_correlation(map_br[:,:,0], image_br[:,:,0])\n",
    "y_offset, x_offset = shift\n",
    "print(y_offset, x_offset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_tr = sat_map[1260-300:1260+300, 15966-300:15966+300, :]\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(cv2.cvtColor(image_tr, cv2.COLOR_BGR2RGB))\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(cv2.cvtColor(map_tr.astype(np.uint8), cv2.COLOR_BGR2RGB))\n",
    "\n",
    "shift, _, _ = phase_cross_correlation(map_br[:,:,0], image_br[:,:,0])\n",
    "y_offset, x_offset = shift\n",
    "print(y_offset, x_offset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_c = sat_map[9688-300:9688+300, 8255-300:8255+300, :]\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(cv2.cvtColor(image_c, cv2.COLOR_BGR2RGB))\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(cv2.cvtColor(map_c.astype(np.uint8), cv2.COLOR_BGR2RGB))\n",
    "\n",
    "shift, _, _ = phase_cross_correlation(map_c[:,:,0], image_c[:,:,0])\n",
    "y_offset, x_offset = shift\n",
    "print(y_offset, x_offset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the transformation between pixel coordinates and UTM coordinates\n",
    "\n",
    "primary = np.array([[900., 900.], #tl\n",
    "                    [492., 18168.], #bl\n",
    "                    [15966., 1260.], #tr\n",
    "                    [15553., 18528.], #br\n",
    "                    [8255., 9688.]]) # c\n",
    "\n",
    "secondary = np.array([[619400., 5736195.],\n",
    "                      [619400., 5734600.],\n",
    "                      [620795., 5736195.],\n",
    "                      [620795., 5734600.],\n",
    "                      [620100., 5735400.]])\n",
    "\n",
    "# Pad the data with ones, so that our transformation can do translations too\n",
    "n = primary.shape[0]\n",
    "pad = lambda x: np.hstack([x, np.ones((x.shape[0], 1))])\n",
    "unpad = lambda x: x[:,:-1]\n",
    "X = pad(primary)\n",
    "Y = pad(secondary)\n",
    "\n",
    "# Solve the least squares problem X * A = Y\n",
    "# to find our transformation matrix A\n",
    "A, res, rank, s = np.linalg.lstsq(X, Y)\n",
    "\n",
    "transform = lambda x: unpad(np.dot(pad(x), A))\n",
    "\n",
    "print(\"Target:\")\n",
    "print(secondary)\n",
    "print(\"Result:\")\n",
    "print(transform(primary))\n",
    "print(\"Max error:\", np.abs(secondary - transform(primary)).max())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}