2_2_pyspatiotemporal_interpol.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[Pyspatiotemporalgeom](http://www.cs.siue.edu/~marmcke/docs/pyspatiotemporalgeom/) based interpolation\n",
    "\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** this library was writen in Python 2. It had to be converted to \n",
    "          Python 3 compatible code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Author: \u001B[1mTiago Ribeiro\u001B[0m\n",
      "\n",
      "Github username: \u001B[1mTiago1Ribeiro\u001B[0m\n",
      "\n",
      "Last updated: 2023-04-04 12:05:14\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.10.9\n",
      "IPython version      : 8.11.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from watermark import watermark\n",
    "print(watermark(author=\"\\033[1m\" + \"Tiago Ribeiro\"+ \"\\033[0m\", \n",
    "                github_username=\"\\033[1m\" + \"Tiago1Ribeiro\"+ \"\\033[0m\", \n",
    "                current_date=True, current_time=True, python=True, \n",
    "                updated=True, iversions=True, globals_= globals())\n",
    "                )"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": "### Setup `pyspatiotemporalgeom`  library for Python 3"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "```bash\n",
    "# Clones the repository with the Python3 version of pyspatiotemporalgeom lib.    \n",
    "git clone https://github.com/CIIC-C-T-Polytechnic-of-Leiria/pyspatiotemporalgeom-for-Python-3.git\n",
    "# Moves pyspatiotemporalgeom folder to the working directory \n",
    "mv pyspatiotemporalgeom-for-Python-3/pyspatiotemporalgeom .\n",
    "# Removes the cloned repository\n",
    "rm -rf pyspatiotemporalgeom-for-Python-3 \n",
    "# on Windows: rmdir /s /q pyspatiotemporalgeom-for-Python-3 \n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:37:12.897142Z",
     "start_time": "2024-05-22T18:37:12.770408Z"
    }
   },
   "source": [
    "import os\n",
    "import re\n",
    "from geomet import wkt\n",
    "from typing import List, Union\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from glob import glob\n",
    "# for debugging purposes\n",
    "from IPython.display import clear_output\n",
    "\n",
    "import sys\n",
    "sys.path.insert(0,\"working_dir\\ESS_Data_Lab_misc\\pyspatiotemporalgeom-for-Python-3\\pyspatiotemporalgeom\")\n",
    "\n",
    "import pyspatiotemporalgeom.structureRegion as structureRegion\n",
    "from pyspatiotemporalgeom.componentMovingRegion import cIntervalRegion\n",
    "from pyspatiotemporalgeom.utilities import hsegLibrary"
   ],
   "outputs": [],
   "execution_count": 3
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data Sources"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:38:29.284514Z",
     "start_time": "2024-05-22T18:38:29.272244Z"
    }
   },
   "source": [
    "DATA_DIR = \"D:\\BurnedAreaUAV_files\\BurnedAreaUAV_dataset\\BurnedAreaUAV_dataset\" # Replace with your data directory# Directory to save interpolated polygons\n",
    "OUT_DIR = \"E://BurnedAreaUAV_files//Interpolation//pyspatiotemporalgeom_interpol\"\n",
    "# WKT containing the manually annotated polygons\n",
    "WKT_FILE = os.path.join(DATA_DIR, 'WKT//train_valid.wkt')\n",
    "# Directory to save PNG format interpolated polygons\n",
    "OUT_DIR_PNG = os.path.join(OUT_DIR, 'PNGs')\n",
    "# Directory to save WKT format interpolated polygons\n",
    "OUT_WKT_FILE = os.path.join(OUT_DIR, \"pyspatial_interpol.wkt\")"
   ],
   "outputs": [],
   "execution_count": 5
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:38:32.680889Z",
     "start_time": "2024-05-22T18:38:31.193403Z"
    }
   },
   "source": [
    "# Convert WKT to list of shapely polygons\n",
    "with open(WKT_FILE, 'r') as wkt_file:\n",
    "    multipolygons_file = wkt_file.readlines()\n",
    "    multipolygons = [wkt.loads(multipolygons_file[i]) for i in range(len(multipolygons_file))]"
   ],
   "outputs": [],
   "execution_count": 6
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Auxiliary functions"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "create_CIR"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:39:15.446341Z",
     "start_time": "2024-05-22T18:39:15.420188Z"
    }
   },
   "source": [
    "def create_CIR(source_SL: List, dest_SL: List, source_time: int, dest_time: int) -> Union[cIntervalRegion, None]:\n",
    "    \"\"\"\n",
    "    Creates a Component Interval Region.\n",
    "    This describes the movement of a moving region over a single time interval.\n",
    "\n",
    "    Args:\n",
    "        source_SL (List): The source segment list.\n",
    "        dest_SL (List): The destination segment list.\n",
    "        source_time (int): The initial instant of the time interval.\n",
    "        dest_time (int): The final instant of the time interval.\n",
    "\n",
    "    Returns:\n",
    "        Union[cIntervalRegion, None]: A Component Interval Region or None.\n",
    "\n",
    "    Note:\n",
    "        source_time and dest_time must be integer timestamps.\n",
    "\n",
    "    Adapted from: https://github.com/most-ieeta/SPT-DataLab/blob/master/pyspatiotemporalgeom/mckenney_final.py\n",
    "    \"\"\"\n",
    "    src_sr = structureRegion.structuralRegion()\n",
    "    src_sr_id = src_sr.addFace(source_SL)\n",
    "\n",
    "    dst_sr = structureRegion.structuralRegion()\n",
    "    dst_sr_id = dst_sr.addFace(dest_SL)\n",
    "\n",
    "    cir = cIntervalRegion()\n",
    "\n",
    "    cir.sourceSR = src_sr\n",
    "    cir.destSR = dst_sr\n",
    "    cir.sourceTime = int(source_time)\n",
    "    cir.destTime = int(dest_time)\n",
    "\n",
    "    if cir.mapComponent(src_sr_id, dst_sr_id):\n",
    "        return cir\n",
    "\n",
    "    return None"
   ],
   "outputs": [],
   "execution_count": 7
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "labelUniqueCycles"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:39:18.529904Z",
     "start_time": "2024-05-22T18:39:18.499624Z"
    }
   },
   "source": [
    "from collections import deque\n",
    "from functools import cmp_to_key\n",
    "def labelUniqueCycles(segs, onlyReturnOuterCycle=False):\n",
    "    '''\n",
    "    will do interior walks around cycles.  \n",
    "\n",
    "    Each cycle will get a unique label.\n",
    "\n",
    "    Label numbers will not necessarily start at 2, and may skip numbers\n",
    "\n",
    "    Will also remove sticks!!!  a full service solution!\n",
    "\n",
    "    Cycles will have unique labels, but holes in particular will have labelling \n",
    "    flipped.  \n",
    "\n",
    "    Use: def switchLabelsForCorrectCycleLabelling( hsegs ): to finalize labels\n",
    "\n",
    "    This function is used to create well-formed regions from possibly \n",
    "    non-wellformed input.\n",
    "\n",
    "    **input**: \n",
    "\n",
    "    segs: a list of segs.\n",
    "\n",
    "    onlyReturnOuterCycle:  if you want a simple region, set this to true!\n",
    "\n",
    "    **output**: a hseg list representing a well formed region. Agian, call def \n",
    "    switchLabelsForCorrectCycleLabelling( hsegs ): to finalize labels!!! \n",
    "     Everything is labeled as an outercycle after this call.\n",
    "    '''\n",
    "    # remove dups from segs\n",
    "    nonDupSegs = []\n",
    "    for s in segs:\n",
    "        if s[0] < s[1]:\n",
    "            nonDupSegs.append(s)\n",
    "        else:\n",
    "            nonDupSegs.append((s[1], s[0]))\n",
    "    seenSegSet = set(nonDupSegs)\n",
    "    segs = list(seenSegSet)\n",
    "    # convert segs to hsegs\n",
    "    hsegs = []\n",
    "    for s in segs:\n",
    "        hsegs.append(((s[0], s[1]), -1, -1))\n",
    "        hsegs.append(((s[1], s[0]), -1, -1))\n",
    "    hsegs.sort(key=cmp_to_key(hsegCompForSorted))\n",
    "    # create a hash table for the seg portions mapped to their index\n",
    "    indexLookup = dict()\n",
    "    for i, h in enumerate(hsegs):\n",
    "        indexLookup[h[0]] = i\n",
    "\n",
    "    # visit unprocessed segs in hseg order.  Each unprocessed will start a new cycle\n",
    "    currLabel = 1\n",
    "    nestedLabel = 2\n",
    "    for i, h in enumerate(hsegs):\n",
    "        if (h[1] == -1 and h[2] == -1) and isLeft(h):  # unprocessed, left\n",
    "            #nestedLabel = currLabel+1\n",
    "            # interior walk the cycle\n",
    "            visitedPoiSet = set()\n",
    "            seenSegSet = set()\n",
    "            visitedHsegStack = deque()\n",
    "            completeVisitedHistoryStack = deque()\n",
    "            currIndex = i\n",
    "            startIndex = i\n",
    "            firstTimeThrough = True\n",
    "            currLabel = nestedLabel\n",
    "            nestedLabel += 1\n",
    "            labellingAbove = True\n",
    "            startHseg = hsegs[startIndex]\n",
    "            while True:\n",
    "                currHseg = hsegs[currIndex]\n",
    "                # check if this is a stick seg.  If so, just label it as invalid and move on\n",
    "                if (currHseg[0][1], currHseg[0][0]) in seenSegSet:\n",
    "                    # we have walked its brother.  this is a stick. label it as such\n",
    "                    hsegs[currIndex] = (currHseg[0], 0, 0)\n",
    "                    brotherIndex = indexLookup[(\n",
    "                        currHseg[0][1], currHseg[0][0])]\n",
    "                    brother = hsegs[brotherIndex]\n",
    "                    brother = (brother[0], 0, 0)\n",
    "                    hsegs[brotherIndex] = brother\n",
    "                    # If this happens to be the brother of the startHseg, we \n",
    "                    # need to unvisit the walked cycle. (since we walked the \n",
    "                    # exterior and may have walked mutliple stick-connected cycles)\n",
    "                    if brother[0] == hsegs[startIndex][0]:\n",
    "                        for itemHseg in completeVisitedHistoryStack:\n",
    "                            if itemHseg[1] != 0:\n",
    "                                itemIndex = indexLookup[itemHseg[0]]\n",
    "                                brotherIndex = indexLookup[(\n",
    "                                    itemHseg[0][1], itemHseg[0][0])]\n",
    "                                hsegs[itemIndex] = (itemHseg[0], -1, -1)\n",
    "                                hsegs[brotherIndex] = (\n",
    "                                    (itemHseg[0][1], itemHseg[0][0]), -1, -1)\n",
    "                        break\n",
    "                # check if we have reached a processed seg (happpnes if we started on a stick)\n",
    "                # we have to process first seg twice\n",
    "                elif (currHseg[1] > 0 or currHseg[2] > 0) and currIndex != startIndex:\n",
    "                    # we must have started on a stick, and followed it to get here.\n",
    "                    # mark everything we followed as such\n",
    "                    visitedHsegStack.appendleft(startHseg)\n",
    "                    for itemHseg in visitedHsegStack:\n",
    "                        itemIndex = indexLookup[itemHseg[0]]\n",
    "                        brotherIndex = indexLookup[(\n",
    "                            itemHseg[0][1], itemHseg[0][0])]\n",
    "                        hsegs[itemIndex] = (itemHseg[0], 0, 0)\n",
    "                        hsegs[brotherIndex] = (\n",
    "                            (itemHseg[0][1], itemHseg[0][0]), 0, 0)\n",
    "                    break\n",
    "                # final case we have an unprocessed, non-stick seg.\n",
    "                # label it appropriately\n",
    "                else:\n",
    "                    # label the current hseg and its brother\n",
    "                    la = currLabel\n",
    "                    lb = -1\n",
    "                    if not labellingAbove:\n",
    "                        lb = currLabel\n",
    "                        la = -1\n",
    "\n",
    "                    hsegs[currIndex] = (currHseg[0], la, lb)\n",
    "                    brotherIndex = indexLookup[(\n",
    "                        currHseg[0][1], currHseg[0][0])]\n",
    "                    brother = hsegs[brotherIndex]\n",
    "                    brother = (brother[0], la, lb)\n",
    "                    hsegs[brotherIndex] = brother\n",
    "\n",
    "                    # if we have closed a loop, pop and fix\n",
    "                    if currHseg[0][0] in visitedPoiSet:\n",
    "                        stopPoi = currHseg[0][0]\n",
    "\n",
    "                        while True:\n",
    "                            # debug statment.  We are about to pop a stack.  if it is empty, this prog\n",
    "                            # will crash.  So, preemptively printout the current region. current seg,\n",
    "\n",
    "                            # end of debug statement\n",
    "                            # pop it, updateLabels\n",
    "                            # print currHseg\n",
    "                            poppedHseg = visitedHsegStack.popleft()\n",
    "\n",
    "                            poppedIndex = indexLookup[poppedHseg[0]]\n",
    "                            poppedIndexBrother = indexLookup[(\n",
    "                                poppedHseg[0][1], poppedHseg[0][0])]\n",
    "                            if hsegs[poppedIndex][1] > 0:\n",
    "                                hsegs[poppedIndex] = (\n",
    "                                    poppedHseg[0], nestedLabel, -1)\n",
    "                                hsegs[poppedIndexBrother] = (\n",
    "                                    (poppedHseg[0][1], poppedHseg[0][0]), nestedLabel, -1)\n",
    "                            elif hsegs[poppedIndex][2] > 0:\n",
    "                                hsegs[poppedIndex] = (\n",
    "                                    poppedHseg[0], -1, nestedLabel)\n",
    "                                hsegs[poppedIndexBrother] = (\n",
    "                                    (poppedHseg[0][1], poppedHseg[0][0]), -1, nestedLabel)\n",
    "                            if poppedHseg[0][0] == stopPoi:\n",
    "                                break\n",
    "                        nestedLabel += 1\n",
    "\n",
    "                # update visited stacks and get the next seg\n",
    "                if firstTimeThrough:\n",
    "                    firstTimeThrough = False\n",
    "                else:\n",
    "                    # if not first time through, record point and seg, check stopping condition\n",
    "                    visitedPoiSet |= set([currHseg[0][0]])\n",
    "                    visitedHsegStack.appendleft(currHseg)\n",
    "                    completeVisitedHistoryStack.appendleft(currHseg)\n",
    "                    if currIndex == startIndex:\n",
    "                        # if we get here, we finished a cycle\n",
    "                        # if we only want the outer cycle, we can jsut return it here\n",
    "                        if onlyReturnOuterCycle:\n",
    "                            hsegs = [h for h in hsegs if h[1] != h[2]]\n",
    "                            return hsegs\n",
    "                        break\n",
    "\n",
    "                seenSegSet |= set([(currHseg[0])])\n",
    "\n",
    "                # get next hseg\n",
    "                prevIndex = currIndex\n",
    "                currIndex = getNextInnerCycleWalkIndex(\n",
    "                    currHseg, hsegs, indexLookup)\n",
    "                # check if we are switching label above\n",
    "                if isLeft(hsegs[prevIndex]) != isLeft(hsegs[currIndex]):\n",
    "                    labellingAbove = not labellingAbove\n",
    "\n",
    "    hsegs = [h for h in hsegs if h[1] != h[2]]\n",
    "    return hsegs"
   ],
   "outputs": [],
   "execution_count": 8
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "get_region_t"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:39:32.230389Z",
     "start_time": "2024-05-22T18:39:32.206116Z"
    }
   },
   "source": [
    "def get_region_t(cir, t):\n",
    "    \"\"\"\n",
    "    Returns the structural region defined at the specified time.\n",
    "\n",
    "    Input: A component interval region and a time.\n",
    "    Output: The structural region defined at the specified time.\n",
    "    \"\"\"\n",
    "\n",
    "    # Extract the structural region defined by cir at time t.\n",
    "    return cir.getStructuralRegionAtTime(t)"
   ],
   "outputs": [],
   "execution_count": 9
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "getStructuralRegionAtTime"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:39:34.293578Z",
     "start_time": "2024-05-22T18:39:34.277998Z"
    }
   },
   "source": [
    "def getStructuralRegionAtTime(self, time):\n",
    "    '''\n",
    "    Extract the structural region defined by this component interval region \n",
    "    at time ``time``.\n",
    "\n",
    "    Call the interpolator, get the triangles, extract the region from those \n",
    "    triangles\n",
    "\n",
    "    **Input**:\n",
    "    time:\n",
    "        The time at which to extract the structural regions\n",
    "\n",
    "    **Returns**:\n",
    "    None:\n",
    "        If the ``time`` is outside the bounds of this interval regions\n",
    "\n",
    "    structural region:\n",
    "        otherwise\n",
    "    '''\n",
    "    # print('1')\n",
    "    if time < self.sourceTime or time > self.destTime:\n",
    "        return None\n",
    "    # print('2')\n",
    "    if time == self.sourceTime:\n",
    "        return self.sourceSR\n",
    "\n",
    "    if time == self.destTime:\n",
    "        return self.destSR\n",
    "    # print('3')\n",
    "    print(\"time ok\", end=\"\\r\")\n",
    "    # find out the percentage of how far along the time interval **time** lies\n",
    "    multiplier = float(time-self.sourceTime) / \\\n",
    "        (self.destTime - self.sourceTime)\n",
    "\n",
    "    # step 1:  get the motion triangles, but maintain the IDs of the\n",
    "    # compnonets involved so we can rebuild the mappings\n",
    "    componentTris = self.getMotionTriangles()\n",
    "    # step 2:  create a structural region to hold the result\n",
    "    sr = structureRegion.structuralRegion()\n",
    "    # print('componentTris', componentTris)\n",
    "    # step 3: go through the triangles for each moving component, extract\n",
    "    # the segs at time **time**, and add the structures\n",
    "    for component in componentTris:\n",
    "        segs = []\n",
    "        for tri in component[2]:\n",
    "            s1 = ((tri[0][0], tri[0][1]), (tri[2][0], tri[2][1]))\n",
    "            s2 = ((tri[1][0], tri[1][1]), (tri[2][0], tri[2][1]))\n",
    "            if tri[2][2] < tri[0][2]:\n",
    "                s1 = (s1[1], s1[0])\n",
    "                s2 = (s2[1], s2[0])\n",
    "            p1x = (multiplier * (s1[1][0] - s1[0][0])) + s1[0][0]\n",
    "            p1y = (multiplier * (s1[1][1] - s1[0][1])) + s1[0][1]\n",
    "            p2x = (multiplier * (s2[1][0] - s2[0][0])) + s2[0][0]\n",
    "            p2y = (multiplier * (s2[1][1] - s2[0][1])) + s2[0][1]\n",
    "            segs.append(((p1x, p1y), (p2x, p2y)))\n",
    "        # now we have the segs.  place them into the return SR\n",
    "        # note, there will never be a line or point in the middle of an interval\n",
    "        if self.sourceSR.getComponentType(component[0]) != 'H':\n",
    "            sr.addFace(segs, component[0])\n",
    "        elif self.sourceSR.getComponentType(component[0]) == 'H':\n",
    "            sr.addHole(segs, self.sourceSR.getIDsForHolesInFace(\n",
    "                component[0]), component[0])\n",
    "        # print('segs: ', segs)\n",
    "    return sr"
   ],
   "outputs": [],
   "execution_count": 10
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "points_list_to_intermediate_wkt"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:39:37.821810Z",
     "start_time": "2024-05-22T18:39:37.796320Z"
    }
   },
   "source": [
    "def points_list_to_intermediate_wkt(points):\n",
    "    \"\"\"\n",
    "\tReturns a string representation of the given list of points in the form:\n",
    "\t\t(x1 y1, ..., xn yn, x1 y1)\n",
    "\t\n",
    "\tInput: A list of points: [(x0, y0), ..., (xn, yn)]\n",
    "\tOutput: (x1 y1, ..., xn yn, x1 y1)\n",
    "\t\n",
    "\tNotes:\n",
    "\tUsed as an intermediate step in the computation of a wkt of a geometry.\n",
    "\n",
    "    Source: https://github.com/most-ieeta/SPT-DataLab/blob/master/pyspatiotemporalgeom/mckenney_final.py\n",
    "    \"\"\"\n",
    "    \n",
    "    if points != None:\n",
    "        _wkt = '('\n",
    "        for point in points:\n",
    "            _wkt += str(point[0]) + ' ' + str(point[1]) + ', '\n",
    "        # Close the cycle.\n",
    "        _wkt += str(points[0][0]) + ' ' + str(points[0][1]) + ')'\n",
    "        return _wkt\n",
    "\n",
    "    return None"
   ],
   "outputs": [],
   "execution_count": 11
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "wkt_to_segment_list"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:39:40.224629Z",
     "start_time": "2024-05-22T18:39:40.203499Z"
    }
   },
   "source": [
    "def wkt_to_segment_list(_wkt):\n",
    "    \"\"\"\n",
    "        Parses a wkt string and returns a list of tuples with line segments \n",
    "        representing a geometry.\n",
    "        \n",
    "        Notes:\n",
    "        Currently the function does not handle holes!\n",
    "        Currently the function handles only POLYGON type geometries!\n",
    "        \n",
    "        Input: A wkt string.\n",
    "        Output: A segment list.\n",
    "        \n",
    "        WKT Examples:\n",
    "        \n",
    "        POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30))\n",
    "        POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))\n",
    "        \n",
    "        MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)), ((15 5, 40 10, 10 20, 5 10, 15 5)))\n",
    "        MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)), ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))\n",
    "        \n",
    "        Source: https://github.com/most-ieeta/SPT-DataLab/blob/master/pyspatiotemporalgeom/mckenney_final.py\n",
    "    \"\"\"\n",
    "    segment_list = []\n",
    "    geom_json = wkt.loads(_wkt)\n",
    "    geometry_type = geom_json['type']\n",
    "    if geometry_type != \"Polygon\":\n",
    "         print('Only Polygon is supported.')\n",
    "         return []\n",
    "\n",
    "    coordinates = geom_json['coordinates']\n",
    "\n",
    "    for i in range(0, len(coordinates)):\n",
    "        values = coordinates[i]\n",
    "        # Face.\n",
    "        if i == 0:\n",
    "            for j in range(0, len(values) - 1):\n",
    "                # Close the segment list.\n",
    "                if (j == len(values) - 2):\n",
    "                    segment_list.append((tuple(values[j]), tuple(values[0])))\n",
    "                else:\n",
    "                    segment_list.append((tuple(values[j]), tuple(values[j+1])))\n",
    "        # Holes.\n",
    "        else:\n",
    "            \"\"\"\n",
    "            for j in range(0, len(values) - 1):\n",
    "                # Close the segment list.\n",
    "                if (j == len(values) - 2):\n",
    "                    segment_list.append((tuple(values[j]), tuple(values[0])))\n",
    "                else:\n",
    "                    segment_list.append((tuple(values[j]), tuple(values[j+1])))\n",
    "            \"\"\"\n",
    "    # Validate the wkt parsing process output.\n",
    "    if segment_list == []:\n",
    "        print(\"Error parsing wkt.\")\n",
    "\n",
    "    return segment_list"
   ],
   "outputs": [],
   "execution_count": 12
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "structural_region_to_wkt"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:39:43.835971Z",
     "start_time": "2024-05-22T18:39:43.819877Z"
    }
   },
   "source": [
    "def structural_region_to_wkt(sregion):\n",
    "    \"\"\"\n",
    "        Gets the wkt of a given structural region.\n",
    "    \"\"\"\n",
    "\n",
    "    # Currently, only POLYGON is being considered.\n",
    "    # In the future, MULTIPOLYGON will also be considered.\n",
    "    geom_wkt = ''\n",
    "    counter = 0\n",
    "\n",
    "    # Iterate through the faces.\n",
    "    for face_id in sregion.F:\n",
    "        face = sregion.F[face_id]\n",
    "\n",
    "        # Labeled segs for this face.\n",
    "        hsegs = hsegLibrary.labelUniqueCycles(face, True)\n",
    "\n",
    "        # List of points, in cyclic order, that define the boundary of the outer cycle of this face.\n",
    "        face_outer_cycle_points = hsegLibrary.getOuterWalkPointSequence(hsegs)\n",
    "\n",
    "        if counter > 0:\n",
    "            geom_wkt += ', '\n",
    "\n",
    "        # Add the intermediate wkt representation of the face with no holes.\n",
    "        geom_wkt += '(' + points_list_to_intermediate_wkt(face_outer_cycle_points)\n",
    "        counter = counter + 1\n",
    "\n",
    "        # Face has holes?\n",
    "        if sregion.F2H != None:\n",
    "            if face_id in sregion.F2H:\n",
    "                #print(\"Has holes.\")\n",
    "                for j in sregion.F2H[face_id]:\n",
    "                    hsegs = hsegLibrary.labelUniqueCycles(sregion.H[j])\n",
    "                    hole_outer_cycle_points = hsegLibrary.getOuterWalkPointSequence(hsegs)\n",
    "                    # Add the intermediate wkt representation of the face's hole(s).\n",
    "                    geom_wkt += ', ' + points_list_to_intermediate_wkt(hole_outer_cycle_points)\n",
    "\n",
    "        # Close the polygon wkt representation.\n",
    "        geom_wkt += ')'\n",
    "\n",
    "    # The geometry is a single polygon (1 face with >= 0 holes).\n",
    "    if counter == 1:\n",
    "        geom_wkt = 'POLYGON ' + geom_wkt\n",
    "    # The geometry has >1 faces each with >= 0 holes.\n",
    "    else:\n",
    "        geom_wkt = 'MULTIPOLYGON (' + geom_wkt + ')'\n",
    "\n",
    "    return geom_wkt"
   ],
   "outputs": [],
   "execution_count": 13
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-22T18:47:05.498216Z",
     "start_time": "2024-05-22T18:47:05.467354Z"
    }
   },
   "source": [
    "def interpolate_polygons(multipolygons_file: List[str], n_samples: int, source_t: int, target_t: int, output_file = \"output.wkt\") -> None:\n",
    "    \"\"\"\n",
    "    Interpolates polygons between consecutive pairs of multipolygons in a file.\n",
    "\n",
    "    Parameters:\n",
    "    - multipolygons_file (List[str]): A list of strings, where each string is a \n",
    "      WKT representation of a multipolygon.\n",
    "    - n_samples (int): The number of intermediate polygons to interpolate \n",
    "      between each pair of consecutive polygons.\n",
    "    - source_t (int): The timestamp of the first multipolygon in the file, \n",
    "      in a format recognized by PySpatioTemporalGeom.\n",
    "    - target_t (int): The timestamp of the last multipolygon in the file, \n",
    "      in a format recognized by PySpatioTemporalGeom.\n",
    "    - output_file (str): The name of the output file to write the interpolated\n",
    "\n",
    "    Example usage:\n",
    "    >>> multipolygons_file = [\"MULTIPOLYGON (((0 0, 0 1, 1 1, 1 0, 0 0)))\", \n",
    "        \"MULTIPOLYGON (((1 1, 1 2, 2 2, 2 1, 1 1)))\"]\n",
    "    >>> interpolated_pol = interpolate_polygons(multipolygons_file, 5, 10, 20)\n",
    "    >>> print(interpolated_pol)\n",
    "    ['POLYGON ((0.0 0.0, 0.0 0.2, 0.2 0.2, 0.2 0.0, 0.0 0.0))', \n",
    "      'POLYGON ((0.0 0.0, 0.0 0.4, 0.4 0.4, 0.4 0.0, 0.0 0.0))', ...] \n",
    "    \"\"\"\n",
    "    polygons_num = len(multipolygons_file)\n",
    "    with open(output_file, 'a') as f:\n",
    "      \n",
    "        # Iterate over all consecutive pairs of polygons\n",
    "        for i in range(polygons_num-1):\n",
    "              f.write(multipolygons_file[0]+ '\\n')\n",
    "              source_geom = wkt_to_segment_list(multipolygons_file[i])\n",
    "              target_geom = wkt_to_segment_list(multipolygons_file[i+1])\n",
    "              # Create the Component Interval Region\n",
    "              cir = create_CIR(source_geom, target_geom, source_t, target_t)\n",
    "              # Time interval and number of samples for this pair of polygons\n",
    "              delta_t = target_t - source_t\n",
    "              dt = delta_t / (n_samples + 1)\n",
    "              \n",
    "              # Interpolate n_samples equally spaced between source and target\n",
    "              for j in range(1, n_samples+1):\n",
    "                  t = j * dt / delta_t + source_t\n",
    "                  # Region at query_time\n",
    "                  region_at = get_region_t(cir, t)\n",
    "                  clear_output()\n",
    "                  region_wkt = structural_region_to_wkt(region_at)\n",
    "                  #f.write(region_wkt + '\\n')\n",
    "        return None"
   ],
   "outputs": [],
   "execution_count": 27
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "invalid_polygons = [180,181,182,193,194,201,202]\n",
    "\n",
    "interpolated_polygons = list()\n",
    "with open(OUT_WKT_FILE, 'w') as file:\n",
    "    file.write('')\n",
    "idx_list = list()\n",
    "\n",
    "for i in range(0, 180):\n",
    "    print(f\"Indice {i}\")\n",
    "    idx_list.append(i)\n",
    "    interpolate_polygons([multipolygons_file[i],multipolygons_file[i+1]], \n",
    "                         n_samples = 99, source_t = 0, target_t = 1, \n",
    "                         output_file = OUT_WKT_FILE)\n",
    "\n",
    "#     99+1+99+1+99 = 399\n",
    "#  |----|----|----|\n",
    "# 179  180  181  182  183\n",
    "interpolate_polygons([multipolygons_file[180],multipolygons_file[183]], \n",
    "                     n_samples = 299, source_t = 0, target_t = 1, \n",
    "                     output_file = OUT_WKT_FILE)\n",
    "                \n",
    "# 18300 \n",
    "\n",
    "for i in range(183, 192):\n",
    "    print(f\"Indice {i}\")\n",
    "    idx_list.append(i)\n",
    "    interpolate_polygons([multipolygons_file[i],multipolygons_file[i+1]], \n",
    "                         n_samples = 99, source_t = 0, target_t = 1, \n",
    "                         output_file = OUT_WKT_FILE)\n",
    "\n",
    "#    99+1+99+1+99 = 299\n",
    "#  |----|----|----|\n",
    "# 192  193  194  195 \n",
    "interpolate_polygons([multipolygons_file[192], multipolygons_file[195]], \n",
    "                     n_samples = 299, source_t = 0, target_t = 1, \n",
    "                     output_file = OUT_WKT_FILE)\n",
    "\n",
    "for i in range(195, 200):\n",
    "    print(f\"Indice {i}\")\n",
    "    idx_list.append(i)\n",
    "    interpolate_polygons([multipolygons_file[i],multipolygons_file[i+1]], \n",
    "                         n_samples = 99, source_t = 0, target_t = 1, \n",
    "                         output_file = OUT_WKT_FILE)\n",
    "\n",
    "# 20000\n",
    "\n",
    "#    99+1+99+1+99 = 299\n",
    "#  |----|----|----|\n",
    "# 200  201  202  203 \n",
    "interpolate_polygons([multipolygons_file[200], multipolygons_file[203]],\n",
    "                     n_samples = 299, source_t = 0, target_t = 1, \n",
    "                     output_file = OUT_WKT_FILE)\n",
    "\n",
    "#20300\n",
    "\n",
    "for i in range(203,  len(multipolygons_file)-1):\n",
    "    print(f\"Indice {i}\")\n",
    "    idx_list.append(i)\n",
    "    interpolate_polygons([multipolygons_file[i],multipolygons_file[i+1]], \n",
    "                         n_samples = 99, source_t = 0, target_t = 1, \n",
    "                         output_file = OUT_WKT_FILE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert WKT to list of shapely polygons\n",
    "with open(OUT_WKT_FILE, 'r') as wkt_file:\n",
    "    interpol_file = wkt_file.readlines()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### wkt2masc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wkt2masc(wkt_file, images_path, orig_dims, out_dims, delete_files=True):\n",
    "    \"\"\" \n",
    "    Converts WKT files to segmentation masks.\n",
    "    Parameters:\n",
    "        wkt_file {str} -- path to the WKT file\n",
    "        images_path {str} -- path to the folder where the masks will be saved\n",
    "        orig_dims {tuple} -- (width, height) original dimensions of the masks \n",
    "        out_dims {tuple} -- (width, height) output dimensions of the masks  \n",
    "    Returns:\n",
    "        Creates PNG images of the masks\n",
    "    \"\"\"\n",
    "\n",
    "    os.makedirs(images_path, exist_ok=True)\n",
    "\n",
    "    if delete_files:\n",
    "        # delete files in the folder, if any\n",
    "        for filename in os.listdir(images_path):\n",
    "            if filename.endswith(\".png\"):\n",
    "                os.remove(os.path.join(images_path, filename))\n",
    "\n",
    "    # open WKT file\n",
    "    wkt = open(wkt_file, 'r')\n",
    "    num_lines = len(wkt.readlines())\n",
    "    cnt = 0\n",
    "    \n",
    "    print(f\"\"\"\n",
    "    {'-'*38}\n",
    "    # \\033[1mProperties of the resulting masks\\033[0m\n",
    "    # Width: {out_dims[0]}, Height: {out_dims[1]}\n",
    "    # Number of masks to create: {num_lines}\n",
    "    {'-'*38}\n",
    "    \"\"\")\n",
    "    \n",
    "    # process each line of the WKT file\n",
    "    with open(wkt_file) as f:\n",
    "        for line in f:\n",
    "            # extract numbers from the line\n",
    "            points = [int(float(number)) for number in re.findall(r'\\d+(?:\\.\\d+)?', line)]\n",
    "            # create empty mask\n",
    "            mask = np.zeros((orig_dims[1],orig_dims[0]), dtype=np.uint8)\n",
    "            # create array with polygon points, with 2 columns (x,y)\n",
    "            arr = np.array(points, dtype=np.int32).reshape((-1,2))\n",
    "            # draw mask\n",
    "            cv2.drawContours(image = mask,\n",
    "                             contours=[arr],\n",
    "                             contourIdx=-1,\n",
    "                             color=(255, 255, 255),\n",
    "                             thickness=-1,  # if > 0, thickness of the contour; if -1, fill object\n",
    "                             lineType=cv2.LINE_AA)\n",
    "            \n",
    "            if out_dims != orig_dims:\n",
    "                # resize frames with Lanczos interpolation\n",
    "                mask = cv2.resize(mask, out_dims, interpolation=cv2.INTER_CUBIC)\n",
    "            # save mask as PNG\n",
    "            cv2.imwrite(os.path.join(images_path, f\"frame_{cnt:06d}.png\"), mask)\n",
    "            cnt += 1\n",
    "            # print progress\n",
    "            print(f\"\\r\\033[1m{cnt}\\033[0m/{num_lines} masks created\", end=\"\\r\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### frames2video"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def frames2video(img_list, nome_ficheiro='video', fps_ = 25, titulo: str = \"\", frame_num_text  = False, font_size: int = 1) -> None:\n",
    "    \"\"\" \n",
    "    Converte lista de imagens em ficheiro AVI com a mesma resolucão da primeira \n",
    "    imagem da lista.\n",
    "      Parametros: - lista de imagens PNG, TIFF, JPEG, BMP, WEBP, STK, LSM ou XCF\n",
    "                  - nome do ficheiro do video\n",
    "      Devolve: salva vídeo no diretório de execucão\n",
    "    \"\"\"\n",
    "    # guarda dimensões da primeira imagem\n",
    "    img = cv2.imread(img_list[0])\n",
    "    height, width, _ = img.shape\n",
    "    size = (width, height)\n",
    "    num_frames =  len(img_list)\n",
    "\n",
    "    img_array = list()\n",
    "    for i in range(len(img_list)):\n",
    "        img = cv2.imread(img_list[i])\n",
    "        img_array.append(img)\n",
    "        print(f\"1. Appending frames {i+1}/{num_frames}\", end=\"\\r\")\n",
    "        \n",
    "    print(\"2. Creating video writer...\", end=\"\\r\")\n",
    "    video = cv2.VideoWriter(filename= nome_ficheiro + '.avi',\n",
    "                            fourcc=cv2.VideoWriter_fourcc(*'mp4v'), fps = fps_,\n",
    "                            frameSize=size)\n",
    "\n",
    "    for i in range(len(img_array)):\n",
    "        if frame_num_text:\n",
    "\n",
    "            frame_number_text = f\"frame_{i:06d}\"\n",
    "            cv2.putText(img_array[i], frame_number_text, (width-300, 50), \n",
    "                            cv2.FONT_HERSHEY_SIMPLEX,font_size, (255, 100, 100), \n",
    "                            2, cv2.LINE_AA)\n",
    "        if titulo:\n",
    "            cv2.putText(img_array[i], titulo, (50, 50), cv2.FONT_HERSHEY_SIMPLEX,\n",
    "                        font_size, (255, 255, 255), 2, cv2.LINE_AA)\n",
    "        \n",
    "        video.write(img_array[i])\n",
    "        print(f\"3. Writing frames to file {i+1}/{num_frames}\", end=\"\\r\")\n",
    "    video.release()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### PNGs and video creation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of lines in the WKT file: 22500\n",
      "\n",
      "    --------------------------------------\n",
      "    # \u001B[1mProperties of the resulting masks\u001B[0m\n",
      "    # Width: 1280, Height: 720\n",
      "    # Number of masks to create: 22500\n",
      "    --------------------------------------\n",
      "    \n",
      "\u001B[1m22500\u001B[0m/22500 masks created\r"
     ]
    }
   ],
   "source": [
    "wkt2masc(wkt_file = OUT_WKT_FILE, images_path = OUT_DIR_PNG, orig_dims = (1280, 720), out_dims = (1280, 720))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3. Writing frames to file 22500/22500\r"
     ]
    }
   ],
   "source": [
    "img_list = sorted(glob(os.path.join(OUT_DIR_PNG, '*.png')))\n",
    "frames2video(img_list = img_list, nome_ficheiro='pyspatial_interpol', \n",
    "             fps_ = 25*10, titulo = \"Pyspatialtemporalgeom interpolation (10x speed)\", \n",
    "             frame_num_text = True, font_size = 1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation of the sampled polygons"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WKT with sampled polygons\n",
    "WKT_FILE_SAMPLED = os.path.join(\"E:/BurnedAreaUAV_files/Interpolation/reference_masks\", \"sampled_masks.txt\")\n",
    "# Directory to save PNG format interpolated polygons\n",
    "OUT_DIR_SAMPLED_PNG = os.path.join(OUT_DIR, 'PNGs_sampled')\n",
    "# create output directory\n",
    "if not os.path.exists(OUT_DIR_SAMPLED_PNG):\n",
    "    os.makedirs(OUT_DIR_SAMPLED_PNG)\n",
    "# Directory to save WKT format interpolated polygons    \n",
    "OUT_WKT_SAMPLED_FILE = os.path.join(OUT_DIR, \"pyspatial_interpol_sampled.wkt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read txt file and extract polygons and their indexes\n",
    "with open(WKT_FILE_SAMPLED, 'r') as f:\n",
    "    polygons = f.readlines()\n",
    "    # extract indexes and polygons\n",
    "    indexes = [int(polygon.split(',')[0]) for polygon in polygons]\n",
    "    polygons_wkt = [polygon.split(',', 1)[1][:-1] for polygon in polygons]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# makes sure the file is empty\n",
    "open(OUT_WKT_SAMPLED_FILE, 'w').close()\n",
    "for i in range(0, len(polygons_wkt)-1):\n",
    "    # calculate number of samples to interpolate\n",
    "    num_samples = (indexes[i+1] - indexes[i])*100 - 1\n",
    "    interpolate_polygons([polygons_wkt[i], polygons_wkt[i+1]], \n",
    "                         n_samples = num_samples, source_t = 0, target_t = 1, \n",
    "                         output_file = OUT_WKT_SAMPLED_FILE)\n",
    "\n",
    "# append last polygon to WKT file\n",
    "with open(OUT_WKT_SAMPLED_FILE, 'a') as f:\n",
    "    f.write(polygons_wkt[-1] + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(OUT_WKT_SAMPLED_FILE, 'a') as f:\n",
    "    f.write(polygons_wkt[-1] + '\\n')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### WKT to PNG conversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    --------------------------------------\n",
      "    # \u001B[1mProperties of the resulting masks\u001B[0m\n",
      "    # Width: 1280, Height: 720\n",
      "    # Number of masks to create: 22501\n",
      "    --------------------------------------\n",
      "    \n",
      "\u001B[1m22501\u001B[0m/22501 masks created\r"
     ]
    }
   ],
   "source": [
    "wkt2masc(wkt_file = OUT_WKT_SAMPLED_FILE, images_path = OUT_DIR_SAMPLED_PNG, orig_dims = (1280, 720), out_dims = (1280, 720))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Video Generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3. Writing frames to file 22501/22501\r"
     ]
    }
   ],
   "source": [
    "img_list = sorted(glob(os.path.join(OUT_DIR_SAMPLED_PNG, '*.png')))\n",
    "frames2video(img_list = img_list, nome_ficheiro='pyspatial_interpol_sampled', \n",
    "             fps_ = 25*10, titulo = \"Pyspatialtemporalgeom interpolation (10x speed)\", \n",
    "             frame_num_text = True, font_size = 1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "interpol",
   "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.9"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}