video_objects_removal.ipynb

{
  "cells":[
    {
      "cell_type":"markdown",
      "source":[
        "### `Minelli Giovanni`\n",
        "# Moving objects removal\n",
        "The project goal is to analyze an input video and return the static background of the scene comprehensive of small dynamic changes over time. The algorithm is supposed to work properly both with input from a static camera and with a free moving camera (dynamic shot). Also, the algorithm works on-the-fly without preassumptions or preprocessing of the video frames."
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "from termcolor import colored\n",
        "import numpy as np\n",
        "import cv2\n",
        "from matplotlib.gridspec import GridSpec\n",
        "from matplotlib import pyplot as plt\n",
        "from scipy import ndimage, misc\n",
        "from skimage.metrics import structural_similarity as ssim\n",
        "from skimage import filters\n",
        "import time\n",
        "import progressbar"
      ],
      "execution_count":1,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "def drawlines(img1,img2,lines,pts1,pts2):\n",
        "    ''' Utility for showing the correspondances between detected points\n",
        "        img1 - image on which we draw the epilines for the points in img2\n",
        "        lines - corresponding epilines '''\n",
        "    r,c = img1.shape\n",
        "    img1 = cv2.cvtColor(img1,cv2.COLOR_GRAY2BGR)\n",
        "    img2 = cv2.cvtColor(img2,cv2.COLOR_GRAY2BGR)\n",
        "    for r,pt1,pt2 in zip(lines,pts1,pts2):\n",
        "        color = tuple(np.random.randint(0,255,3).tolist())\n",
        "        x0,y0 = map(int, [0, -r[2]\/r[1] ])\n",
        "        x1,y1 = map(int, [c, -(r[2]+r[0]*c)\/r[1] ])\n",
        "        # img1 = cv2.line(img1, (x0,y0), (x1,y1), color,1)\n",
        "        img1 = cv2.circle(img1,tuple(pt1),5,color,-1)\n",
        "        img2 = cv2.circle(img2,tuple(pt2),5,color,-1)\n",
        "    return img1,img2"
      ],
      "execution_count":2,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "<img src=\"docs\/auto_canny.png\" height=\"400\" \/>"
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "def auto_canny(image, sigma=0.33):\n",
        "\t''' Apply Canny edge detector and thresholding using a lower and upper boundary on the gradient values.\n",
        "\tLower value of sigma indicates a tighter threshold, whereas a larger value of sigma gives a wider threshold '''\n",
        "\n",
        "\t# compute the median of the single channel pixel intensities\n",
        "\tv = np.median(image)\n",
        "\t# apply automatic Canny edge detection using the computed median\n",
        "\tlower = int(max(0, (1.0 - sigma) * v))\n",
        "\tupper = int(min(255, (1.0 + sigma) * v))\n",
        "\tedged = cv2.Canny(image, lower, upper)\n",
        "\t# return the edged image\n",
        "\treturn edged"
      ],
      "execution_count":3,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "- https:\/\/docs.opencv.org\/4.5.1\/d5\/d0f\/tutorial_py_gradients.html\n",
        "- https:\/\/en.wikipedia.org\/wiki\/Image_gradient\n",
        "\n",
        "`\"The gradient at each image point is a 2D vector with the components given by the derivatives in the horizontal and vertical directions. Each gradient vector points in the direction of largest possible intensity increase.\"`\n",
        "<img src=\"docs\/count_oriented.png\" height=\"800\" \/>\n",
        "\n",
        "An high threshold is applied on the given image to keep only the point were an high gradient is found and therefore we suspect the presence of the border of the object. The input image is an intermediary step of the `spot_the_diff` function"
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "def count_oriented(img, show=False):\n",
        "  ''' Count the points in a binary image taking in consideration the orientation '''\n",
        "\n",
        "    # sobel derivatives\n",
        "  derivX = cv2.Sobel(img, cv2.CV_32F, 1, 0)\n",
        "  derivY = cv2.Sobel(img, cv2.CV_32F, 0, 1)\n",
        "\n",
        "  mag = cv2.magnitude(np.absolute(derivX), np.absolute(derivY)) # absolute since a limitation of opencv from docs\n",
        "  \n",
        "    # thresholding of the magnitude values\n",
        "  thresh = 1000\n",
        "  _, mask = cv2.threshold(mag, thresh, 255, cv2.THRESH_BINARY)\n",
        "  mask = np.uint8(mask>0)\n",
        "  \n",
        "  if show:\n",
        "    plt.figure(figsize=(15,10)); plt.title(\"Mask magnitude min: {} max: {}\".format(np.int(np.min(mag[mag>0])), np.int(np.max(mag))))\n",
        "    plt.imshow(cv2.cvtColor(mask*255, cv2.COLOR_BGR2RGB))\n",
        "    plt.show()\n",
        "  \n",
        "  return np.sum(mask)"
      ],
      "execution_count":4,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "def expand_borders(m): return cv2.dilate(m, None, iterations=5)\n",
        "def enhance_borders(im, k=np.array([[-1,-1,-1], [-1, 9,-1], [-1,-1,-1]])): return cv2.GaussianBlur(im, (3, 3), 0) # cv2.filter2D(im, -1, k) #"
      ],
      "execution_count":5,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "**This functions is used to determine if the bound (the next_bound) around the moving object is done correctly evaluating the difference between the two images in the ipotetical static zone.**\n",
        "<img src=\"docs\/distances.png\" height=\"1000\" \/>\n",
        "Confronting different functions of distance and configuration for thresholding or enhancing the differences i decided to use the structural similarity for a low points detection and for the aim of disambiguate over a possible equality of the images. In case of an evident presence of a difference between the two a more \"precise\" function is used to get the shape of the moving object.\n",
        "- absdiff (double check the equality with an higher thresh value)\n",
        "- Morphological transformation: close \n",
        "- Morphological transformation: gradient on the border of the shape found\n",
        "- Canny of both images using the gradient as mask\n",
        "- count_oriented of edges found by Canny and sum\n",
        "\n",
        "The image with the higher value is supposed to contain \"the object\". The return value is unbalanced by 0.01% to prioritize the False negative and minimize the False positive.\n",
        "\n",
        "<img src=\"docs\/spot_the_diff.png\" height=\"1000\" \/>"
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "index = 0\n",
        "\n",
        "def spot_the_diff(image1_gray, image2_gray, log=False, show_work=True, debug_init=True):\n",
        "    ''' Given two images return the one with the object  \n",
        "          image1_gray: [prev_bound - next_bound] * prev_show\n",
        "          image2_gray: [prev_bound - next_bound] * next_frame\n",
        "    ret=> if bound was correct the diff function will find dirt in image1 (=>1) or image1 and image2 are equals (=>0)\n",
        "          if bound was incorrect the diff function will contain some moving parts of image2 (=>2)\n",
        "             (can happen that also identifies small part like light changes of image1 but returns 2) ''' \n",
        "\n",
        "    if debug_init:\n",
        "        global index\n",
        "        plt.figure(figsize=(15,15))\n",
        "        plt.subplot(121); plt.title(\"PREV {}\".format(index))\n",
        "        plt.imshow(cv2.cvtColor(image1_gray, cv2.COLOR_BGR2RGB))\n",
        "        plt.subplot(122); plt.title(\"NEXT\")\n",
        "        plt.imshow(cv2.cvtColor(image2_gray, cv2.COLOR_BGR2RGB))\n",
        "        plt.show()\n",
        "        plt.imsave(\"spot_debug\/prev_{}.jpg\".format(index),cv2.cvtColor(image1_gray, cv2.COLOR_BGR2RGB))\n",
        "        plt.imsave(\"spot_debug\/next_{}.jpg\".format(index),cv2.cvtColor(image2_gray, cv2.COLOR_BGR2RGB))\n",
        "\n",
        "    (score, sim) = ssim(image1_gray, image2_gray, full=True)\n",
        "    sim = np.uint8(sim*255)\n",
        "    thresh_sim =  np.uint8(np.logical_not(cv2.threshold(sim, 205, 255, cv2.THRESH_BINARY)[1])*255)\n",
        "    thresh_sim = cv2.medianBlur(thresh_sim, 7)\n",
        "    sum_sim = np.sum(thresh_sim>0)\n",
        "\n",
        "      # lower bound (low thresh, high detection)\n",
        "    if sum_sim<100:\n",
        "        print(colored(\"Equality_l\", 'blue'))\n",
        "        return 0\n",
        "\n",
        "    dif = cv2.absdiff(image1_gray, image2_gray)\n",
        "    thresh_dif = cv2.medianBlur(dif, 5)\n",
        "\n",
        "      # higher bound (high thresh, low detection)\n",
        "    if not np.any(thresh_dif>50):\n",
        "        print(colored(\"Equality_h\", 'blue'))\n",
        "        return 0\n",
        "\n",
        "    if show_work:\n",
        "        plt.figure(figsize=(20,20))\n",
        "        plt.subplot(131); plt.title(\"ssim\")\n",
        "        plt.imshow(cv2.cvtColor(np.uint8(thresh_sim), cv2.COLOR_BGR2RGB))\n",
        "\n",
        "    adaptive=0\n",
        "    detected = 110\n",
        "    a=None\n",
        "    while(detected>80 if detected<100 else detected>105):\n",
        "        a = np.uint8(thresh_dif>adaptive)*255\n",
        "        detected = ((np.sum(a>0)\/sum_sim)*100).round(3)\n",
        "        if log: print(\"Not {}%\".format(detected.round(2)))\n",
        "        adaptive+=10\n",
        "\n",
        "    mask = cv2.morphologyEx(a,cv2.MORPH_CLOSE, np.ones((12,12),np.uint8))\n",
        "    mask = cv2.medianBlur(mask,3)\n",
        "    gradient = cv2.morphologyEx(mask, cv2.MORPH_GRADIENT, np.ones((6,6),np.uint8))\n",
        "\n",
        "    if show_work:\n",
        "        plt.subplot(132); plt.title(\"Mask closed [{}]\".format(adaptive-10))\n",
        "        plt.imshow(cv2.cvtColor(np.uint8(mask), cv2.COLOR_BGR2RGB))\n",
        "        plt.subplot(133); plt.title(\"Mask gradient [{}%]\".format(detected.round(1)))\n",
        "        plt.imshow(cv2.cvtColor(np.uint8(gradient), cv2.COLOR_BGR2RGB))\n",
        "        plt.show()\n",
        "\n",
        "      # canny in the border of the movement areas\n",
        "    blurred1 = enhance_borders(image1_gray*(gradient>0))\n",
        "    blurred2 = enhance_borders(image2_gray*(gradient>0))\n",
        "\n",
        "    edges1 = auto_canny(blurred1)\n",
        "    edges2 = auto_canny(blurred2)\n",
        "\n",
        "      # count of canny oriented pixels in the border of the movement areas \n",
        "    i1=count_oriented(edges1)\n",
        "    i2=count_oriented(edges2)\n",
        "    if log: print(i1)\n",
        "    if log: print(i2)\n",
        "\n",
        "    if show_work:\n",
        "        plt.figure(figsize=(20,20))\n",
        "        plt.subplot(121)\n",
        "        plt.title(\"Canny 1 [{}]\".format(i1))\n",
        "        plt.imshow(cv2.cvtColor(edges1, cv2.COLOR_BGR2RGB))\n",
        "        plt.subplot(122)\n",
        "        plt.title(\"Canny 2 [{}]\".format(i2))\n",
        "        plt.imshow(cv2.cvtColor(edges2, cv2.COLOR_BGR2RGB))\n",
        "        plt.show()\n",
        "\n",
        "    if show_work:\n",
        "        plt.figure(figsize=(20,20))\n",
        "        plt.subplot(121); plt.title(\"index {} sim_score {}\".format(index, (score*100).round(1), 15), fontsize=20); \n",
        "        plt.imshow(cv2.cvtColor(image1_gray*(mask>0), cv2.COLOR_BGR2RGB))\n",
        "        plt.subplot(122)\n",
        "        plt.imshow(cv2.cvtColor(image2_gray*(mask>0), cv2.COLOR_BGR2RGB))\n",
        "\n",
        "      # prioritize the false negative cases to minimize the false positives \n",
        "    if i1>(i2+(0.01*i2)):\n",
        "        if show_work: plt.title(\"Oggetto in 1\", fontsize=20); plt.show()\n",
        "        print (colored(\"#Oggetto in 1 [{}]\".format(i1-i2),'blue')); return 1;\n",
        "    else:\n",
        "        if show_work: plt.title(\"Oggetto in 2\", fontsize=20); plt.show()\n",
        "        print(colored(\"#Oggetto in 2 [{}]\".format(i2-i1),'blue')); return 2;"
      ],
      "execution_count":31,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "**Detection of the camera movement (prev frame to next frame) by SIFT feature points.**\n",
        "\n",
        "Assuming a static camera the points found in the background image will always be the same. If there is a good amount of cordinates matching between the background image and the next frame then no waping will be perfomed.\n",
        "<img src=\"docs\/sift_moving.png\" height=\"600\" \/>\n",
        "\n",
        "\n",
        "Otherways recalculate foundamental matrix from eight-point method using only background feature points (computing the distance of point from epipolar line <a href=\"https:\/\/www.mdpi.com\/2076-3417\/10\/1\/268\">[Jung, S.; Cho, Y.; Kim, D.; Chang, M. Moving Object Detection from Moving Camera...]<\/a>)\n",
        "<img src=\"docs\/point_line_distance_formula.png\" \/>\n",
        "Moreover can be applied a RANSAC algorithm to clean result from outliers. And then compute the image transformation of both the previous frame and the background image."
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "des_recovery = None\n",
        "\n",
        "def warp(prev_show, prev_frame, frame, show_work=False, show_result=False):\n",
        "      # find the keypoints and descriptors with SIFT\n",
        "    global des_recovery\n",
        "    global mean_prev\n",
        "\n",
        "    if not des_recovery is None:\n",
        "        kp1, des1 = des_recovery\n",
        "    else:\n",
        "        kp1, des1 = sift.detectAndCompute(prev_show,None)\n",
        "        des_recovery = (kp1, des1)\n",
        "    kp2, des2 = sift.detectAndCompute(frame,None)\n",
        "    \n",
        "      # determine if camera is static, \n",
        "    if len(kp1)!=0 and len(kp2)!=0:\n",
        "      pp1 = np.array([[int(p.pt[0]),int(p.pt[1])] for p in kp1])\n",
        "      pp2 = np.array([[int(p.pt[0]),int(p.pt[1])] for p in kp2])\n",
        "      matches = int(np.sum((pp1[:,None] == pp2).all(2).any(1))\/len(kp1)*100)\n",
        "      if(matches<THRESH_STATIC_OR_MOVING): # the prev images need warping and old SIFT points are trashed\n",
        "          print( colored(\"Warping\", \"cyan\") )\n",
        "          des_recovery = None\n",
        "      else: # no need to warp\n",
        "          return (prev_show, prev_frame, frame)\n",
        "\n",
        "    matches1 = []; matches1_bg = []\n",
        "    matches2 = []; matches2_bg = []\n",
        "      # match the descriptors and build the corrspondences arrays of good matches\n",
        "    matches_pairs = flann.knnMatch(des1, des2,k=2)\n",
        "    for i,(m,n) in enumerate(matches_pairs):\n",
        "        if m.distance < 0.7*n.distance:\n",
        "            matches1.append(kp1[m.queryIdx].pt)\n",
        "            matches2.append(kp2[m.trainIdx].pt)\n",
        "    matches1 = np.asarray(matches1, np.int32)\n",
        "    matches2 = np.asarray(matches2, np.int32)\n",
        "\n",
        "    if len(matches1)>MIN_MATCH_COUNT:\n",
        "        F, mask = cv2.findFundamentalMat(matches1,matches2,cv2.FM_LMEDS)\n",
        "          # we select only inlier points\n",
        "        matches1 = matches1[mask.ravel()==1]\n",
        "        matches2 = matches2[mask.ravel()==1]\n",
        "\n",
        "          # find epilines corresponding to points in left image (first image) and drawing its lines on right image\n",
        "        lines2 = cv2.computeCorrespondEpilines(matches1.reshape(-1,1,2),1,F)\n",
        "        lines2 = lines2.reshape(-1,3)\n",
        "\n",
        "        if show_work:\n",
        "            img3,img4 = drawlines(prev_show,frame,lines2,matches2,matches1)\n",
        "            plt.figure(figsize=(20,10))\n",
        "            plt.subplot(121); plt.title(\"lines {0} points {1}\".format(len(lines2), len(matches1))); plt.imshow(img3)\n",
        "            plt.subplot(122); plt.imshow(img4)\n",
        "            plt.show()\n",
        "\n",
        "          # keep only points near <1px to an epipolar line\n",
        "        for pt1,pt2 in zip(matches1,matches2):\n",
        "            pixels_err = PIXELS_ERR\n",
        "            for line in lines2:\n",
        "                pixels_err = min(pixels_err,abs(line[0] * pt1[0] + line[1] * pt1[1] + line[2]) \/ np.sqrt((line[0] * line[0]) + (line[1]*line[1])))\n",
        "            if pixels_err<1:\n",
        "                matches1_bg.append(pt1)\n",
        "                matches2_bg.append(pt2)\n",
        "\n",
        "        print(\"bg matches: \",len(matches1_bg),\"\/\",len(matches1))\n",
        "\n",
        "    if len(matches1_bg)>MIN_MATCH_COUNT:\n",
        "          # retrieving perspective homography\n",
        "        src = np.float32(matches1_bg)\n",
        "        dst = np.float32(matches2_bg)\n",
        "        H, masked = cv2.findHomography(src,dst, cv2.RANSAC, 5.0)\n",
        "        dst_show = cv2.warpPerspective(np.uint16(prev_show)+1,H,(prev_frame.shape[1], prev_frame.shape[0]))\n",
        "        dst = cv2.warpPerspective(np.uint16(prev_frame)+1,H,(prev_frame.shape[1], prev_frame.shape[0]))\n",
        "        \n",
        "          # correct the last pixel at margins\n",
        "        dst_show[cv2.dilate(np.uint8(dst_show==0), None)>0]=0\n",
        "        dst[cv2.dilate(np.uint8(dst==0), None)>0]=0\n",
        "\n",
        "          # warp\n",
        "        prev_show = np.uint8(((dst_show-1)*(dst_show>0))+((dst_show<1)*frame))\n",
        "        prev_frame = np.uint8(((dst-1)*(dst>0))+((dst<1)*frame))\n",
        "        mean_prev = [prev_show]\n",
        "\n",
        "        if show_result:\n",
        "            plt.figure(figsize=(20,10))\n",
        "            plt.subplot(121); plt.title.set_text('Left warp')\n",
        "            plt.imshow(cv2.cvtColor(dst,cv2.COLOR_GRAY2RGB))\n",
        "            plt.subplot(122); plt.title.set_text('Compelete prev frame warped')\n",
        "            plt.imshow(cv2.cvtColor(prev_frame,cv2.COLOR_GRAY2RGB))\n",
        "            plt.show()\n",
        "        \n",
        "        return (prev_show, prev_frame, frame)\n",
        "    else:\n",
        "        print( colored(\"Not enough matches are found - {}\/{}\".format(len(matches1_bg), MIN_MATCH_COUNT), \"red\") )\n",
        "        mean_prev = [prev_show]\n",
        "\n",
        "        return None"
      ],
      "execution_count":21,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "**FLOW**"
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "The basic understanding of the algorithms is the idea of recognize objects in movement confronting just two frames at each cycle. Found a moving object with a proper diff strategy we can confront the bound of the oject in movement with the bound found at the cycle before and then make assumptions on the movement: (↓↑←→): in a direction, (B): backwards or (F): forward.\n",
        ".  \n",
        "- (↓↑←→) In this case we can gain \"good pixels\" both from the previous frame (in the B zone) and from the next frame (in the A zone)\n",
        "\n",
        "<img src=\"docs\/sliding.png\" \/>\n",
        "\n",
        "- (B) In this case we can take pixels from the next frame in the zone of previous movement (A zone) \n",
        "\n",
        "<img src=\"docs\/backward.png\" \/>\n",
        "\n",
        "- (C) all the movement zone remain unused\n",
        "\n",
        "Those regions of improvement can suffer of bad bounds detection and also of the presence of old pixels of moving objects in the background image to show (due to initialization with the the first frame). To disambiguate between those situations the function `spot_the_diff` is used to confront the prev_show (the background image) and the next_frame. The response can be \"equality\" (under a thresh value) and so we can use the newer pixels for the update, or it can indicate the image in which the object is.\n",
        "\n",
        "Moreover all the other static pixels (not detected as movement in the two frames difference) are used for the update of the background.\n"
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "''' Hyperp to tune properly '''\n",
        " # FLANN parameters\n",
        "FLANN_INDEX_KDTREE = 1\n",
        "index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5)\n",
        "search_params = dict(checks=50)\n",
        "sift = cv2.SIFT_create()\n",
        "flann = cv2.FlannBasedMatcher(index_params,search_params)\n",
        " # WARP parameters\n",
        "THRESH_STATIC_OR_MOVING = 12 # %\n",
        "PIXELS_ERR = 1\n",
        "MIN_MATCH_COUNT = 10\n",
        " #other parameters\n",
        "sim_thresh = 165\n",
        "min_rect_body_ratio_thresh = 0.005\n",
        "min_rect_thresh = (5, 10)\n",
        "max_bounds_per_zone = 10\n",
        "frames_to_skip = 2\n",
        "mean_prev_size = 5\n",
        "\n",
        "show_work = False\n",
        "debug_spot = False\n",
        "show_result = False\n",
        "print_log = True\n",
        "\n",
        " # iterators\n",
        "index = 0\n",
        "prev_bounds = []\n",
        "frame = None\n",
        "prev_frame = None\n",
        "des_recovery = None\n",
        "output_video = []\n",
        "mean_prev = []\n",
        "        \n",
        "\n",
        "# video_name=files[np.random.randint(0,len(files))]\n",
        "video_name=\"ex\/test.avi\"\n",
        "\n",
        "cap = cv2.VideoCapture(video_name)\n",
        "tot_frames = 0\n",
        "while True:\n",
        "    for _ in range(frames_to_skip): ret, frame_BGR = cap.read() # skip\n",
        "    ret, frame_BGR = cap.read()\n",
        "    if frame_BGR is None: break\n",
        "    tot_frames+=1\n",
        "cap.release()\n",
        "\n",
        "cap = cv2.VideoCapture(video_name)\n",
        "bar = progressbar.ProgressBar(maxval=tot_frames, widgets=[progressbar.Bar('=', '[', ']'), ' ', progressbar.Percentage()])\n",
        "bar.start()\n",
        "# for i in range(100):\n",
        "#     for _ in range(frames_to_skip): ret, frame_BGR = cap.read() # skip\n",
        "#     ret, frame_BGR = cap.read()\n",
        "while True:\n",
        "    for _ in range(frames_to_skip): ret, frame_BGR = cap.read() # skip\n",
        "    ret, frame_BGR = cap.read()\n",
        "    if frame_BGR is None: break\n",
        "    \n",
        "    # restart\n",
        "    start = time.time()\n",
        "    index+=1\n",
        "    prev_frame = frame\n",
        "    frame = cv2.cvtColor(np.uint8(frame_BGR),cv2.COLOR_BGR2GRAY)\n",
        "    if index<2:\n",
        "        min_rect_body_ratio_thresh = (frame.shape[0] * frame.shape[1] ) * min_rect_body_ratio_thresh\n",
        "        prev_show = frame.copy()\n",
        "        mean_prev.append(prev_show.copy())\n",
        "        continue\n",
        "    handle_camera_movement = warp(prev_show, prev_frame, frame)\n",
        "    if handle_camera_movement is None:\n",
        "        frame = prev_frame\n",
        "        continue\n",
        "\n",
        "    prev_show, prev_frame, next_frame = handle_camera_movement\n",
        "    \n",
        "      # identify movement areas in prev_frame to next_frame \n",
        "    (score, sim) = ssim(prev_frame, cv2.GaussianBlur(next_frame, (3, 3), 0), full=True)\n",
        "    sim = np.uint8(sim*255)\n",
        "    thresh = cv2.threshold(sim, sim_thresh, 1, cv2.THRESH_BINARY)[1]\n",
        "    thresh = np.uint8(np.logical_not(thresh))\n",
        "    thresh = cv2.medianBlur(thresh, 11) # remove small pieces\n",
        "    thresh = cv2.morphologyEx(thresh,cv2.MORPH_CLOSE, np.ones((9,9),np.uint8)) # unify pieces\n",
        "\n",
        "    static_area = np.uint8(np.logical_not(thresh))\n",
        "    movement_area = thresh\n",
        "\n",
        "      # get rid of the frame margins\n",
        "    if np.any(movement_area[0,10:-10]) and not np.any(movement_area[10,10:-10]):\n",
        "        movement_area[:10,:] = 0\n",
        "    if np.any(movement_area[-1,10:-10]) and not np.any(movement_area[-10,10:-10]):\n",
        "        movement_area[-10:-1,:] = 0\n",
        "    if np.any(movement_area[10:-10,0]) and not np.any(movement_area[10:-10,10]):\n",
        "        movement_area[:,:10] = 0\n",
        "    if np.any(movement_area[10:-10,-1]) and not np.any(movement_area[10:-10,-10]):\n",
        "        movement_area[:,-10:-1] = 0\n",
        "    \n",
        "      # crete rectangles from countours\n",
        "    rects = []\n",
        "    bounds = []\n",
        "    recovery_bounds = []\n",
        "    contours, hierarchy = cv2.findContours(movement_area,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)\n",
        "    bounded_image=movement_area*prev_frame\n",
        "    for c in contours:\n",
        "        x,y,w,h = cv2.boundingRect(c)\n",
        "        if w > min_rect_thresh[0] and h > min_rect_thresh[1]:\n",
        "            m=np.zeros_like(movement_area, np.uint8)\n",
        "            m[y:y+h,x:x+w] = 1\n",
        "            rects.append({\"x\":x,\"y\":y,\"w\":w,\"h\":h,\"m\":m})\n",
        "        \n",
        "    if print_log: print(\"rects: \", len(rects))\n",
        "\n",
        "      # create indipendent bounds from rectangles\n",
        "    while len(rects)!=0:\n",
        "        x = rects[0][\"x\"];y = rects[0][\"y\"];w = rects[0][\"w\"];h = rects[0][\"h\"]\n",
        "        m = rects[0][\"m\"]\n",
        "        m_expanded = expand_borders(m) \n",
        "        cv2.rectangle(bounded_image,(x,y),(x+w,y+h),(255,255,0),2)\n",
        "        rects.remove(rects[0])\n",
        "\n",
        "        shape_modified = True\n",
        "        while shape_modified:\n",
        "            shape_modified = False\n",
        "            for r_in in rects:\n",
        "                if np.any(np.logical_and(m_expanded, r_in[\"m\"])):\n",
        "                    rects.remove(r_in)\n",
        "                    cv2.rectangle(bounded_image,(r_in[\"x\"],r_in[\"y\"]),(r_in[\"x\"]+r_in[\"w\"],r_in[\"y\"]+r_in[\"h\"]),(115,115,0),1)\n",
        "                    m[r_in[\"y\"]:r_in[\"y\"]+r_in[\"h\"],r_in[\"x\"]:r_in[\"x\"]+r_in[\"w\"]] = 1\n",
        "                    shape_modified = True\n",
        "                    break\n",
        "        bounds.append(m)\n",
        "\n",
        "    if print_log: print(\"bounds: \", len(bounds))\n",
        "\n",
        "      # get rid of the distorsion due to low movement or stopping object\n",
        "    for prev_bound in prev_bounds:\n",
        "        movement_zones = []\n",
        "        for idx_bound in range(len(bounds)):\n",
        "            if np.array_equal(np.logical_and(prev_bound, bounds[idx_bound]), prev_bound):\n",
        "                movement_zones.append(idx_bound)\n",
        "        if len(movement_zones)>=max_bounds_per_zone: # (ignore it) push forward the previous bound to the next round\n",
        "            for idx_zone in reversed(movement_zones):\n",
        "                bounds.pop(idx_zone)\n",
        "            recovery_bounds.append(prev_bound)\n",
        "        elif len(movement_zones)>=2: # (use it) group and add the rectangle\n",
        "            canvas = np.zeros(movement_area.shape)\n",
        "            \n",
        "            for idx_bound in reversed(movement_zones):\n",
        "                canvas[bounds[idx_bound]>0] = 1\n",
        "                bounds.pop(idx_bound)\n",
        "\n",
        "            (y, x) = np.where(canvas == 1)\n",
        "            (topy, topx) = (np.min(y), np.min(x))\n",
        "            (bottomy, bottomx) = (np.max(y), np.max(x))\n",
        "            cv2.rectangle(bounded_image,(topx,topy),(bottomx,bottomy),(115,115,0),3)\n",
        "            canvas[topx:bottomx+1, topy:bottomy+1] = 1\n",
        "            bounds.append(canvas)\n",
        "\n",
        "      # improvement areas \n",
        "    a = [] # areas to ignore into zone [prev_bounds - next_bounds]. The remaining is updated with next_frame\n",
        "    b = [] # areas to update into zone [next_bounds - prev_bounds]. Update with prev_frame\n",
        "    recover_prev_bound = [0]*len(prev_bounds)\n",
        "    if not (prev_bounds is None) or len(prev_bounds) != 0:\n",
        "        for bound in bounds:\n",
        "            a_current = []\n",
        "            b_current = []\n",
        "            recover_current_bound = []\n",
        "            for prev_bound_id,prev_bound in enumerate(prev_bounds):\n",
        "                if (np.sum(prev_bound[bound==1])\/np.sum(bound))>0.5 or (np.sum(bound[prev_bound==1])\/np.sum(prev_bound))>0.5:\n",
        "                    if print_log: print(\"--cycle-- bound matching\")\n",
        "                    if np.any(np.logical_and(prev_bound, bound)):\n",
        "                        if print_log: print(\"crossing\")\n",
        "\n",
        "                        static_prev_zone = np.uint8(np.logical_and(prev_bound, np.logical_not(bound)))\n",
        "                        moving_next_zone = np.logical_and(np.logical_not(prev_bound), bound)\n",
        "                        # if not np.all(prev_bound[bound==1]) and not np.all(bound[prev_bound==1]):\n",
        "                        if np.any(static_prev_zone) and np.any(moving_next_zone):\n",
        "                            imm = spot_the_diff(static_prev_zone*prev_show, static_prev_zone*next_frame, log=print_log, show_work=show_work, debug_init=debug_spot)\n",
        "                            if print_log: print(\"is moving\")\n",
        "                            if imm == 0 or imm == 1: # (update) images equality or obj residual in prev_show \n",
        "                                a_current.append(bound)\n",
        "                                b.append(np.uint8(np.logical_not(moving_next_zone)))\n",
        "                                recover_current_bound.append(bound)\n",
        "                            else: # (ignore it) probably the bound wasn't recognized properly\n",
        "                                area_to_ignore = np.uint8(np.logical_or(prev_bound,bound))\n",
        "                                a_current.append(area_to_ignore)\n",
        "                                recover_current_bound.append(area_to_ignore)\n",
        "                        elif np.any(static_prev_zone):\n",
        "                            imm = spot_the_diff(static_prev_zone*prev_show, static_prev_zone*next_frame, log=print_log, show_work=show_work, debug_init=debug_spot)\n",
        "                            if print_log: print(\"is getting smaller\")\n",
        "                            if imm == 0 or imm == 1: # (update) \n",
        "                                if print_log: print(\"and there is good stuff\")\n",
        "                                a_current.append(bound)\n",
        "                                recover_current_bound.append(bound)\n",
        "                            else:\n",
        "                                a_current.append(prev_bound) # (ignore it)\n",
        "                                recover_current_bound.append(prev_bound)\n",
        "                        else:\n",
        "                            if print_log: print(\"is getting bigger\")\n",
        "                            a_current.append(bound) # (ignore it)\n",
        "                            recover_current_bound.append(bound)\n",
        "            if len(a_current)>0:\n",
        "                a_all = np.all(a_current, axis=0)\n",
        "                a.append(a_all)\n",
        "                if len(b_current)>0:\n",
        "                    b_all = np.all(b_current, axis=0)\n",
        "                    b.append(np.logical_and(np.logical_not(a_all), b_all))\n",
        "            if len(recover_current_bound)>0:\n",
        "                rec_all = np.all(recover_current_bound, axis=0) \n",
        "                contours, hierarchy = cv2.findContours(np.uint8(rec_all),cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)\n",
        "                if len(contours)>1:\n",
        "                    for idx in range(len(contours)):\n",
        "                        canvas = np.zeros_like(next_frame, np.uint8)\n",
        "                        cv2.drawContours(canvas, contours, idx, 1, -1)\n",
        "                        recovery_bounds.append(canvas)\n",
        "                else:\n",
        "                    recovery_bounds.append(rec_all)\n",
        "            \n",
        "    if show_result:\n",
        "        plt.figure(figsize=(20,10))\n",
        "        plt.subplot(131); plt.title(\"Static area estim\")\n",
        "        plt.imshow(cv2.cvtColor(np.uint8(static_area*next_frame), cv2.COLOR_GRAY2RGB))\n",
        "        plt.subplot(132); plt.title(\"Movement area estim\")\n",
        "        plt.imshow(cv2.cvtColor(np.uint8(movement_area*prev_frame), cv2.COLOR_GRAY2RGB))\n",
        "\n",
        "        # apply good static pixels\n",
        "    for old_zone in b:\n",
        "        prev_show[old_zone>0]=(old_zone*prev_show)[old_zone>0]\n",
        "    for non_static_zone in a:\n",
        "        static_area[non_static_zone>0] = 0\n",
        "    if len(a)>0: prev_show[static_area>0]=(static_area*next_frame)[static_area>0]\n",
        "\n",
        "    if show_result:\n",
        "        plt.subplot(132); plt.title(\"Bounds\")\n",
        "        plt.imshow(cv2.cvtColor(np.uint8(bounded_image), cv2.COLOR_GRAY2RGB))\n",
        "        plt.title(\"Show\")\n",
        "        plt.subplot(133); plt.imshow(cv2.cvtColor(np.uint8(prev_show), cv2.COLOR_GRAY2RGB))\n",
        "        plt.show()\n",
        "        \n",
        "    mean_prev.append(prev_show.copy())\n",
        "    if len(mean_prev)>mean_prev_size:\n",
        "        mean_prev.pop(0)\n",
        "\n",
        "    output_video.append(np.hstack([frame.copy(),np.mean(mean_prev, axis=0).copy() if mean_prev_size>0 else prev_show.copy()]))\n",
        "    \n",
        "    #log\n",
        "    print(colored(\"bounds: {0} - prev_bounds: {1} - recovery_bounds: {2}\".format(len(bounds),len(prev_bounds),len(recovery_bounds)),\"yellow\"))\n",
        "    print(colored(\"time: {}\".format(time.time()-start),\"yellow\"))\n",
        "    print()\n",
        "    \n",
        "    #restart\n",
        "    prev_bounds = recovery_bounds if len(prev_bounds) > 0 else bounds # ERR\n",
        "    bar.update((index-1)+1)\n",
        "\n",
        "cap.release()\n",
        "bar.finish()\n",
        "out = cv2.VideoWriter('project.avi',cv2.VideoWriter_fourcc(*'DIVX'), 30, (frame.shape[1]*2,frame.shape[0]))\n",
        "\n",
        "for i in range(len(output_video)):\n",
        "    out.write(cv2.cvtColor(np.uint8(output_video[i]), cv2.COLOR_GRAY2RGB))\n",
        "    out.write(cv2.cvtColor(np.uint8(output_video[i]), cv2.COLOR_GRAY2RGB))\n",
        "out.release()"
      ],
      "execution_count":30,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    }
  ],
  "metadata":{
    
  },
  "nbformat":4,
  "nbformat_minor":0
}