18110062_3D_CV_without_depth_warping.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "id": "D680Bt01VS0N"
   },
   "outputs": [],
   "source": [
    "# Importing necessary Libraries\n",
    "import cv2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [],
   "source": [
    "#____________ Choose the Images Dataset numer for stitching the images _____________#\n",
    "#######################\n",
    "Image_Dataset_No = \"000003345\"\n",
    "# <----- from [\"000000029\", \"000000292\", \"000000705\", \"000001524\", \"000002812\", \"000003345\"]\n",
    "#######################\n",
    "\n",
    "#____________ Choose if want to use in-built functions or not  _____________#\n",
    "#######################\n",
    "in_built = False       # <-----\n",
    "#######################\n",
    "\n",
    "flag = \"\"\n",
    "if(in_built==True):\n",
    "    flag = \"_in_built\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper Functions:\n",
    "#### A) Detect and extract key features from a given image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [],
   "source": [
    "def detect_and_extract_features(I):\n",
    "    #_____ Using SIFT Algorithm from the in-built cv2 library\n",
    "    sift_algo = cv2.xfeatures2d.SIFT_create()\n",
    "    features, f_descriptors = sift_algo.detectAndCompute(I, None)\n",
    "    return features, f_descriptors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B) Feature Matching [Simple Brute Force Method] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [],
   "source": [
    "def match_features(f1, f2, d1, d2, I1, I2):\n",
    "    #_____ Using Brute Force Matching Algo from the in-built cv2 library that matches the features from the given two images\n",
    "    matching_algo = cv2.BFMatcher(cv2.NORM_L2, True)\n",
    "    matchings = matching_algo.match(d1, d2)\n",
    "#     matchings_sort = sorted(matchings, key = lambda x:x.distance)\n",
    "#     match_I = cv2.drawMatches(I1,f1,I2,f2,matchings_sort[:50],outImg = None, flags=2)\n",
    "#     cv2.imwrite('Matches_I1_I2.png', match_I)\n",
    "    return matchings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C) Homography Matrix Calculation using Best Point Correspondences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_homography(point_correspondences):\n",
    "    #___ Creating a transformation matrix that will stitch the given 3 images based on their point_correspondences\n",
    "    # Implementation of the Direct Linear Transformation\n",
    "    A = [] # Assemble\n",
    "    for points in point_correspondences:\n",
    "        x1, y1, x2, y2 = points\n",
    "        A.append([x1, y1, 1, 0, 0, 0, -x2*x1, -x2*y1, -x2])\n",
    "        A.append([0, 0, 0, x1, y1, 1, -y2*x1, -y2*y1, -y2])        \n",
    "    A = np.asarray(A)\n",
    "    U, E, Vt = np.linalg.svd(A) # Using Singular Value Decomposition to find the tranformation vector [Last column of Vt]\n",
    "    H = (Vt[-1,:] / Vt[-1,-1]).reshape(3, 3) # Normalization\n",
    "    return H # 3*3 Homography Transformation Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### D) Using RANSAC (RANdom SAmple Consensus) Algorithm that best fits the Feature matchings of two images\n",
    "- Focuses largely on the Inliners"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ransac_algo(point_correspondences, threshold, N=500):\n",
    "    maxInliers = [] # Stores the maximum Inliers\n",
    "    Homography_matrix = None\n",
    "    \n",
    "    # Run this for multiple iterations to find the best line fit that has maximum no. of inliers\n",
    "    for i in range(N):\n",
    "        # Randomly select 4 points to estimate a fit line\n",
    "        indices = np.random.choice(range(len(point_correspondences)),4, replace = False)\n",
    "        points = map(point_correspondences.__getitem__, indices) \n",
    "        H = find_homography(points)\n",
    "        inliers = [] # To calculate the no. of inliners using this fit\n",
    "        for points in point_correspondences:\n",
    "            x1, y1, x2, y2 = points\n",
    "            vector1 = np.asarray([x1, y1, 1]).T\n",
    "            transformed_vector = H.dot(vector1) # Find the tranformed vector\n",
    "            transformed_vector /= transformed_vector[-1] #Normalization\n",
    "            vector2 = np.asarray([x2,y2,1])\n",
    "            # Inlier if the point lies very close inside the threshold area.\n",
    "            if(np.linalg.norm(transformed_vector-vector2)<=threshold):\n",
    "                inliers.append(points)\n",
    "        if(len(inliers) > len(maxInliers)):\n",
    "            maxInliers = inliers\n",
    "            Homography_matrix = H\n",
    "    # return the best estimated homography matrix\n",
    "    return Homography_matrix, maxInliers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main Function:\n",
    "#### Process:\n",
    "    - Detect and Extract Features from both the Images\n",
    "    - Feature Matching and Homography Matrix estimation using RANSAC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(I1,I2, in_built):\n",
    "    f1,d1 = detect_and_extract_features(I1) # For image-1\n",
    "#     I1_dash = cv2.drawKeypoints(I1, f1,outImage = None)\n",
    "#     cv2.imwrite('sift_I1.png', I1_dash)\n",
    "    print(\"No. of features extracted using SIFT from Ia: {}\".format(len(f1)))\n",
    "    f2,d2 = detect_and_extract_features(I2) # For image-2\n",
    "#     I2_dash = cv2.drawKeypoints(I2, f2,outImage = None)\n",
    "#     cv2.imwrite('sift_I2.png', I2_dash)\n",
    "    print(\"No. of features extracted using SIFT from Ib: {}\".format(len(f2)))\n",
    "\n",
    "    matchings = match_features(f1, f2, d1, d2, I1, I2)\n",
    "    print(\"No. of features matched using individual SIFT descriptors: {}\".format(len(matchings)))\n",
    "    point_correspondences = []\n",
    "    pcAs = []\n",
    "    pcBs = []\n",
    "    for match in matchings:\n",
    "        x1, y1 = f1[match.queryIdx].pt\n",
    "        x2, y2 = f2[match.trainIdx].pt\n",
    "        pcAs.append(f1[match.queryIdx].pt)\n",
    "        pcBs.append(f2[match.trainIdx].pt)\n",
    "        point_correspondences.append([x1, y1, x2, y2])\n",
    "\n",
    "    threshold = 5\n",
    "    \n",
    "    ### Without the use of inbuilt functions\n",
    "    if(in_built == False):\n",
    "        H, maxx = ransac_algo(point_correspondences,threshold,800)\n",
    "        print(\"Max Inliers\", len(maxx))\n",
    "    else:### Homography estimation using inbuilt functions\n",
    "        H, status = cv2.findHomography( np.matrix(pcAs), np.matrix(pcBs), cv2.RANSAC, threshold)\n",
    "\n",
    "#     print(\"Homography Matrix:\", H)\n",
    "    \n",
    "    return H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just for Matching_Example\n",
    "# I2 = cv2.imread(\"./2_2.jpg\")\n",
    "# I1 = cv2.imread(\"./2_3.jpg\")\n",
    "# H_12 = main(I1,I2, in_built)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading the Datasets and Calculating all the requires Homography matricies for Stitching 4 images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No. of features extracted using SIFT from Ia: 1720\n",
      "No. of features extracted using SIFT from Ib: 1723\n",
      "No. of features matched using individual SIFT descriptors: 759\n",
      "Max Inliers 376\n",
      "No. of features extracted using SIFT from Ia: 1723\n",
      "No. of features extracted using SIFT from Ib: 1605\n",
      "No. of features matched using individual SIFT descriptors: 704\n",
      "Max Inliers 276\n",
      "No. of features extracted using SIFT from Ia: 1605\n",
      "No. of features extracted using SIFT from Ib: 1652\n",
      "No. of features matched using individual SIFT descriptors: 699\n",
      "Max Inliers 317\n"
     ]
    }
   ],
   "source": [
    "#Loading Datasets\n",
    "\n",
    "I4 = cv2.imread(f\"./RGBD dataset/{Image_Dataset_No}/im_0.jpg\")\n",
    "I3 = cv2.imread(f\"./RGBD dataset/{Image_Dataset_No}/im_1.jpg\")\n",
    "I2 = cv2.imread(f\"./RGBD dataset/{Image_Dataset_No}/im_2.jpg\")\n",
    "I1 = cv2.imread(f\"./RGBD dataset/{Image_Dataset_No}/im_3.jpg\")\n",
    "\n",
    "# The  main function here does the feature detection, extraction and returns a generated homography matrix\n",
    "# For stitching 4 images we use the following homographies to stitch the images and build the panaroma!\n",
    "H_12 = main(I1,I2, in_built) # Homography of I1 with respect to I2\n",
    "H_23 = main(I2,I3, in_built) # Homography of I2 with respect to I3\n",
    "H_34 = main(I3,I4, in_built) # Homography of I3 with respect to I4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WARPING:\n",
    "#### Process:\n",
    "- Create a template for the final stitched image\n",
    "- Set the Base image into the template [Using the 2nd image of the Panaroma as the base]\n",
    "- Generate the warp of the left-most image\n",
    "- Stitch the warped left-most image to the template, meanwhile using the Pyramid blending for better results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Template frame Sizes\n",
    "h_max = (I1.shape[0] + I2.shape[0] + I3.shape[0] + I4.shape[0])\n",
    "w_max = (I1.shape[1] + I2.shape[1] + I3.shape[1] + I4.shape[1])*2\n",
    "stitched_image = np.zeros((h_max,w_max,3))\n",
    "# Offset for the Base Image\n",
    "h_offset = h_max//2 - I3.shape[0]//2\n",
    "w_offset = w_max//3 - I3.shape[1]//2\n",
    "## Base Image set in the to be stitched_image\n",
    "stitched_image[h_offset : h_offset+I3.shape[0], w_offset : w_offset+I3.shape[1], :] = I3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████| 443/443 [00:04<00:00, 100.52it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Stitching I4 to I3: [Leftmost Image to the base image]\n",
    "\n",
    "temp_image = np.zeros((h_max,w_max,3)) # To store the Warped I4 Image\n",
    "common_area_xs = [] # To find the width and position of the intersecting region\n",
    "\n",
    "x_min = np.min(np.where(stitched_image != 0)[1]) # Minimum location of the already stitched image\n",
    "\n",
    "#_________ The commented method below is to find the Warped image I4 that can be stitched to I3 easily.\n",
    "#_________ But this method isn't effective because it finds the transformed I4 using the H34 homography and does not\n",
    "#_________ necessarily have integral pixel locations. Thus blank dots found in the output Image.\n",
    "\n",
    "################-----------------------##################\n",
    "# for y in range(I4.shape[0]):\n",
    "#     for x in range(I4.shape[1]):\n",
    "#         v_dash = np.linalg.inv(H_34).dot([x,y,1])\n",
    "#         v_dash /= v_dash[2]\n",
    "#         x_dash, y_dash = int(v_dash[0]), int(v_dash[1])\n",
    "        \n",
    "#         temp_image[y_dash+h_offset-1 : y_dash+h_offset+2, x_dash+w_offset-1 : x_dash+w_offset+2,:] = I4[y,x,:]\n",
    "        \n",
    "#         if(x_dash+w_offset > x_min):\n",
    "#             common_area_xs.append(x_dash)\n",
    "################-----------------------##################\n",
    "\n",
    "#### Alternate method is to first calculate the expected bounding box of the Warped I4 image and inverse map each pixel \n",
    "#___ of the bounding box to an tranformed pixel in I4 image. This ensures regularity in the Output Image.\n",
    "\n",
    "# Bounding box limits calculation\n",
    "bounds_I4 = [ [0,0,1], [0,I4.shape[0]-1,1], [I4.shape[1]-1,0,1], [I4.shape[1]-1,I4.shape[0]-1,1] ]\n",
    "Xs, Ys = [], []\n",
    "H_inv = np.linalg.inv(H_34)\n",
    "for point in bounds_I4:\n",
    "    v_dash = H_inv.dot(point) # Warping of the boundary points\n",
    "    v_dash /= v_dash[2]\n",
    "    x_dash, y_dash = int(v_dash[0]), int(v_dash[1])\n",
    "    Xs.append(x_dash)\n",
    "    Ys.append(y_dash)\n",
    "    \n",
    "# Find reverse mapping for all pixels lying in Min to Max limits:\n",
    "width_limits = range(max(min(Xs),-w_offset), min(max(Xs)+1,w_max-w_offset))\n",
    "height_limits = range(max(min(Ys), -h_offset), min(max(Ys)+1,h_max-h_offset))\n",
    "\n",
    "for y in tqdm(height_limits):\n",
    "    for x in width_limits:\n",
    "        v_dash = H_34.dot([x, y, 1]) # Reverse map to Image I4\n",
    "        v_dash = v_dash/v_dash[2]\n",
    "        x_dash, y_dash = int(v_dash[0]), int(v_dash[1])\n",
    "        \n",
    "        # Important Condition to check: the inverse-mapped pixel should lie in I4:\n",
    "        if(x_dash >= 0 and x_dash < I4.shape[1] and y_dash >= 0 and y_dash < I4.shape[0] ): \n",
    "            stitched_image[y + h_offset, x + w_offset,:] = I4[y_dash,x_dash,:]\n",
    "        if(x + w_offset >= x_min): # Calculating the Intersection Area for Blending\n",
    "            common_area_xs.append(x)\n",
    "cv2.imwrite(\"./results/{}/stitch_no_depth{}_1_2.jpg\".format(Image_Dataset_No,flag),stitched_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████| 418/418 [00:04<00:00, 100.99it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Stitching I2 to I3: [3rd Right Image to the base image]\n",
    "\n",
    "temp_image = np.zeros((h_max,w_max,3)) # To store the Warped I2 Image\n",
    "common_area_xs = [] # To find the width and position of the intersecting region\n",
    "\n",
    "x_max = np.max(np.where(stitched_image != 0)[1]) # Maximum width location of the already stitched image\n",
    "\n",
    "#_________ The commented method below is to find the Warped image I2 that can be stitched to I3 easily.\n",
    "#_________ But this method isn't effective because it finds the transformed I2 using the H23 homography and does not\n",
    "#_________ necessarily have integral pixel locations. Thus blank dots found in the output Image.\n",
    "\n",
    "################-----------------------##################\n",
    "# for y in range(I2.shape[0]):\n",
    "#     for x in range(I2.shape[1]):\n",
    "#         v_dash = H_23.dot([x,y,1])\n",
    "#         v_dash /= v_dash[2]\n",
    "#         x_dash, y_dash = int(v_dash[0]), int(v_dash[1])      \n",
    "#         temp_image[y_dash+h_offset-1 : y_dash+h_offset+2, x_dash+w_offset-1 : x_dash+w_offset+2,:] = I2[y,x,:]\n",
    "################-----------------------##################\n",
    "\n",
    "#### Alternate method is to first calculate the expected bounding box of the Warped I2 image and inverse map each pixel \n",
    "#___ of the bounding box to an tranformed pixel in I2 image. This ensures regularity in the Output Image.\n",
    "\n",
    "# Bounding box limits calculation\n",
    "bounds_I2 = [ [0,0,1], [0,I2.shape[0]-1,1], [I2.shape[1]-1,0,1], [I2.shape[1]-1,I2.shape[0]-1,1] ]\n",
    "Xs, Ys = [], []\n",
    "for point in bounds_I2:\n",
    "    v_dash = H_23.dot(point) # Warping of the boundary points\n",
    "    v_dash /= v_dash[2]\n",
    "    x_dash, y_dash = int(v_dash[0]), int(v_dash[1])\n",
    "    Xs.append(x_dash)\n",
    "    Ys.append(y_dash)\n",
    "# Find reverse mapping for all pixels lying in Min to Max limits:\n",
    "width_limits = range(max(min(Xs),-w_offset), min(max(Xs)+1,w_max-w_offset))\n",
    "height_limits = range(max(min(Ys), -h_offset), min(max(Ys)+1,h_max-h_offset))\n",
    "H_inv = np.linalg.inv(H_23)\n",
    "for y in tqdm(height_limits):\n",
    "    for x in width_limits:\n",
    "        v_dash = H_inv.dot([x, y, 1]) # Reverse map to Image I2\n",
    "        v_dash = v_dash/v_dash[2]\n",
    "        x_dash, y_dash = int(v_dash[0]), int(v_dash[1])\n",
    "        # Important Condition to check: the inverse-mapped pixel should lie in I2:\n",
    "        if(x_dash >= 0 and x_dash < I2.shape[1] and y_dash >= 0 and y_dash < I2.shape[0] ):\n",
    "            stitched_image[y + h_offset, x + w_offset,:] = I2[y_dash,x_dash,:]\n",
    "        if(0 <= x + w_offset < x_max):\n",
    "            common_area_xs.append(x) # Calculating the Intersection Area for Blending\n",
    "cv2.imwrite(\"./results/{}/stitch_no_depth{}_1_2_3.jpg\".format(Image_Dataset_No,flag),stitched_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████████████████████████████████████████████| 488/488 [00:07<00:00, 66.51it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Stitching I1 to I3: [Farthest image to the base image]\n",
    "\n",
    "temp_image = np.zeros((h_max,w_max,3)) # To store the Warped I1 Image\n",
    "common_area_xs = [] # To find the width and position of the intersecting region\n",
    "\n",
    "x_max = np.max(np.where(stitched_image != 0)[1]) # Maximum width location of the already stitched image\n",
    "\n",
    "#_________ The commented method below is to find the Warped image I2 that can be stitched to I3 easily.\n",
    "#_________ But this method isn't effective because it finds the transformed I2 using the H23 homography and does not\n",
    "#_________ necessarily have integral pixel locations. Thus blank dots found in the output Image.\n",
    "\n",
    "################-----------------------##################\n",
    "# for y in range(I1.shape[0]):\n",
    "#     for x in range(I1.shape[1]):\n",
    "#         v_dash = ((H_12).dot(H_23)).dot([x,y,1])\n",
    "#         v_dash /= v_dash[2]\n",
    "#         x_dash, y_dash = int(v_dash[0]), int(v_dash[1])         \n",
    "#         temp_image[y_dash+h_offset-2 : y_dash+h_offset+3, x_dash+w_offset-2 : x_dash+w_offset+3,:] = I1[y,x,:]\n",
    "################-----------------------##################\n",
    "\n",
    "#### Alternate method is to first calculate the expected bounding box of the Warped I1 image and inverse map each pixel \n",
    "#___ of the bounding box to an tranformed pixel in I1 image. This ensures regularity in the Output Image.\n",
    "\n",
    "# Bounding box limits calculation\n",
    "bounds_I1 = [ [0,0,1], [0,I1.shape[0]-1,1], [I1.shape[1]-1,0,1], [I1.shape[1]-1,I1.shape[0]-1,1] ]\n",
    "Xs, Ys = [], []\n",
    "HH = (H_12).dot(H_23)\n",
    "for point in bounds_I1:\n",
    "    v_dash = HH.dot(point) # Warping of the boundary points\n",
    "    v_dash = v_dash/v_dash[2]\n",
    "    x_dash, y_dash = int(v_dash[0]), int(v_dash[1])\n",
    "    Xs.append(x_dash)\n",
    "    Ys.append(y_dash)\n",
    "# Find reverse mapping for all pixels lying in Min to Max limits:\n",
    "width_limits = range(max(min(Xs),-w_offset), min(max(Xs)+1,w_max-w_offset))\n",
    "height_limits = range(max(min(Ys), -h_offset), min(max(Ys)+1,h_max-h_offset))\n",
    "\n",
    "H_inv = np.linalg.inv((H_12).dot(H_23))\n",
    "for y in tqdm(height_limits):\n",
    "    for x in width_limits:\n",
    "        v_dash = H_inv.dot([x, y, 1]) # Reverse map to Image I1\n",
    "        v_dash = v_dash/v_dash[2]\n",
    "        x_dash, y_dash = int(v_dash[0]), int(v_dash[1])\n",
    "        # Important Condition to check: the inverse-mapped pixel should lie in I1:\n",
    "        if(x_dash >= 0 and x_dash < I1.shape[1] and y_dash >= 0 and y_dash < I1.shape[0] and x + w_offset < w_max and y + h_offset < h_max):\n",
    "            stitched_image[y + h_offset, x + w_offset,:] = I1[y_dash,x_dash,:]\n",
    "        if(0 <= x + w_offset < x_max):\n",
    "            common_area_xs.append(x)\n",
    "cv2.imwrite(\"./results/{}/Panorama_no_depth{}_1_2_3_4.jpg\".format(Image_Dataset_No,flag),stitched_image)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "18110062_3D-CV_Assignment_1.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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": 1
}